A Historical-Ecological Approach to Understanding the Effects of Timber Harvest on
Wa’uumst (Devil’s Club, Oplopanax horridus) in Lax’yip Madii Lii

by
Adrian Smith
BPl Natural Resource Planning, UNBC, 2019

THESIS SUBMITTED IN PARTIAL FULFILLMENT OF
THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE
IN
NATURAL RESOURCES AND ENVIRONMENTAL STUDIES

UNIVERSITY OF NORTHERN BRITISH COLUMBIA
August 2023
© Adrian Smith 2023

Abstract
Industrial land-use has had profound impacts on Indigenous peoples’ homelands throughout
Canada. Over the last century, logging practices in British Columbia have severed peoples’
connections to the land-base, creating access challenges and disrupting the availability of
important plant resources. Devil’s club (Oplopanax horridus) is an important medicinal plant to
communities in northwestern British Columbia, and there are mounting concerns about the
impacts of logging on this ethnobotanically salient species.
Over the last 70 years, swathes of productive forests throughout Gitxsan homelands have
been impacted by the logging industry. Wilp (house) Luutkudziiwus has seen ~12% of their
territory of Madii Lii altered by industrial-scale clearcut logging, which has left a mosaic of
even-aged cutblocks that have potentially altered devil’s club habitat. To investigate and detail
the potential impacts of logging on devil’s club in Madii Lii territory, this research measured
devil’s club health and vigour across a chrono-sequence of clearcut-logged sites compared with
traditionally managed and unlogged sites in the territory. In addition, Species Distribution
Models in conjunction with two different climate change scenarios (SSP245 moderate GHGs and
SSP585 high GHGs) were used to predict habitat suitability for devil’s club in Madii Lii and
throughout British Columbia to the end of the 21st century.
Results suggest that, compared to traditionally managed areas (and controls), extensive
logging in the territory negatively impacted the health and vigour of devil’s club, especially the
most desirable individuals (with the largest and thickest stems). In addition, the predictive
modelling for both climate scenarios suggest that climate change will increasingly impact highquality devil’s club habitat.

ii

Table of Contents
Abstract ......................................................................................................................................... ii
Table of Contents ......................................................................................................................... iii
List of Tables .................................................................................................................................. v
List of Figures ...............................................................................................................................vi
Acknowledgement / Dedication ................................................................................................. vii
Chapter 1 ........................................................................................................................................ 1
1.0. INTRODUCTION .............................................................................................................. 1
1.1. Background ......................................................................................................................... 6
1.1.1. Ethnobotanical Context ................................................................................................. 6
1.1.2. Logging Context .......................................................................................................... 11
1.2. Research Objectives ......................................................................................................... 14
1.3. Species Distribution Models Literature Review ............................................................ 17
Chapter 2 ...................................................................................................................................... 21
2.0. RESEARCH METHODS ................................................................................................ 21
2.1. Introduction ...................................................................................................................... 21
2.2. Study Area ......................................................................................................................... 22
2.3. Site Selection ..................................................................................................................... 25
2.4. Data Collection.................................................................................................................. 29
2.4.1. Ecological Variables .................................................................................................... 30
2.4.2. Presence-Absence Data ............................................................................................... 31
2.4.3. Devil’s Club Measurements ........................................................................................ 31
2.5. Data Analysis .................................................................................................................... 32
2.6. Species Distribution Modelling ....................................................................................... 33
2.6.1. Madii Lii Predicted Habitat Model .............................................................................. 36
2.6.2. British Columbia Predictive Habitat Model ................................................................ 38
Chapter 3 ...................................................................................................................................... 40
3.0. RESULTS .......................................................................................................................... 40
3.1. Devil’s Club Morphometric Analysis ............................................................................. 40
3.1.1. Abundance, Stem Length, Stem Diameter .................................................................. 40

iii

3.1.2. Examine Relationship between Shade and Leaf Area on Stem Length ...................... 44
3.1.3. Examine Relationship between Shade and Leaf Area on Stem Diameter ................... 46
3.1.4. Relationships between Devil’s Club Metrics .............................................................. 47
3.1.5. Relationship between Elevation and Devil’s Club Metrics ......................................... 48
3.2. Canopy Openness and Leaf Area Analysis .................................................................... 52
3.2.1. Effect of Shade and leaf area on abundance ................................................................ 53
3.3. Tree Cover analysis .......................................................................................................... 55
3.4. Species Distribution Modelling Results .......................................................................... 56
3.4.1. Response plot............................................................................................................... 57
3.4.2. Interpretation of response plots ................................................................................... 59
3.4.3. BIOCLIM Model Evaluation ...................................................................................... 59
3.4.4. General Linear Regression method for the SDM ........................................................ 60
3.4.5. Average Score Model .................................................................................................. 61
3.5. Madii lii Predictive Habitat Suitability Generation ...................................................... 63
3.5.1. Baseline BIOCLIM SDM profile ................................................................................ 63
3.5.2. MaxEnt Model ............................................................................................................. 64
3.5.3. Random Forest Model ................................................................................................. 67
3.6. Predictive Habitat Suitability for Devil’s Club in Madii Lii 2021–2100 ..................... 68
3.6.1. Training and Testing datapoints with BIOCLIM Model ............................................. 71
3.6.2. Maxent Suitability Model ............................................................................................ 72
3.6.3. Random Forest and SVM Model ................................................................................. 74
3.6.4. Random Forest Predicted Habitat Suitability Model 2021–2100................................ 76
Chapter 4 ...................................................................................................................................... 79
4.0. DISCUSSION ........................................................................................................................ 79
4.1. Looking back ..................................................................................................................... 79
4.2. Looking Forward .............................................................................................................. 86
4.3. Limitations of Research ................................................................................................... 92
Chapter 5 ...................................................................................................................................... 94
5.0. CONCLUSION ..................................................................................................................... 94
5.1 Recommendations .............................................................................................................. 95
5.2 Future Research................................................................................................................. 97
Literature Cited ........................................................................................................................... 99
Appendix .................................................................................................................................... 123

iv

List of Tables
Table 1. Thesis research objectives and the research questions and objectives ............................ 16
Table 2. Ecological criteria, parameters and to determine the location of control sampling sites 26
Table 3. Pairwise results from Dunn post-hoc tests ...................................................................... 42
Table 4. Summary of BIOCLIM model evaluation....................................................................... 60
Table 5. Summary of model evaluation for the BIOCLIM test-train model. ................................ 64
Table 6. Summary of model evaluation scores for MaxEnt model ............................................... 67
Table 7. Summary of model evaluation scores for Random Forest model ................................... 67
Table 8. Summary of model evaluation scores for BIOCLIM model evaluation. ........................ 71
Table 9. Summary of MaxEnt model evaluation........................................................................... 73
Table 10. Summary of Random Forest model evaluation. ............................................................ 75

v

List of Figures
Figure 1. Devil’s club green/white flower, budding in early June. ................................................. 7
Figure 2. Devil's club red berries ripening in mid-August. ............................................................. 7
Figure 3. Madii Lii and the Northwest coast of British Columbia. ............................................... 22
Figure 4. Luutkudziiwus’Xsi Gwin Hauums and Madii Lii territories. ........................................ 23
Figure 5. Devil's club sampling sites on Luutkudziiwus traditional territory Madii Lii. .............. 28
Figure 6. Devil’s club species occurrence within Madii Lii ........................................................ 34
Figure 7. British Columbia recorded locations of devil’s club. .................................................... 37
Figure 8. Distribution of median stem abundance per 100 m2 (0.01 ha) ....................................... 41
Figure 9. Distribution of median stem length of devil’s club plants in different stand ages......... 43
Figure 10. Distribution of stem diameter of devil’s club plants different stand ages.................... 44
Figure 11. Relationship between devil’s club stem length and canopy closure. ........................... 45
Figure 12. Relationship between devil’s club stem length and devil’s club leaf area................... 45
Figure 13. Regression relationship between devil’s club stem diameter and canopy shade. ........ 46
Figure 14. Relationship between devil’s club stem diameter and devil’s club leaf area............... 47
Figure 15. Relationship of devil’s club stem abundance and elevation ........................................ 48
Figure 16. Relationship of devil’s club stem length (cm) and elevation ....................................... 49
Figure 17. Relationship of devil's club stem diameter (mm) and elevation .................................. 50
Figure 18. Regression analysis for devil’s club stem length and stand age . ................................ 51
Figure 19. Regression analysis for devil's club stem abundance and stand age. ........................... 51
Figure 20. Relationship of devil's club stem diameter and stand age . .......................................... 52
Figure 21. Distribution of the median leaf area (cm2) of devil’s club plants stand ages. ............. 53
Figure 22. Relationship of devil’s club stem abundance and available shade . ............................ 54
Figure 23. Relationship between devil’s club stem abundance and devil’s club leaf area............ 54
Figure 24. Frequency of occurrences of tree species where devils club was present.................... 55
Figure 25. Training and testing points generated for BIOCLIM SDM. ........................................ 56
Figure 26. Multiple response plots for devil’s club probability of occurrence. ............................ 58
Figure 27. SDM habitat suitability values from BIOCLIM model and presence/absence model. 60
Figure 28. SDM continuous habitat suitability values from GLM model and presence/absence
model. ............................................................................................................................................ 61
Figure 29. SDM using BIOCLIM, Domain, GLM and SVM model profiles ............................... 62
Figure 30. Madii Lii devil’s club Average Score SDM BIOCLIM, Domain, GLM and SVM .... 62
Figure 31. Madii Lii PHM for BIOCLIM model for the SSP126 2021–2040. ............................ 63
Figure 32. Response plots for the environmental variables for MaxEnt model generation. ......... 65
Figure 33. PHM Madii Lii for the SSP245:2021–2040 MaxEnt layer.......................................... 66
Figure 34. PHM Madii Lii for SSP245:2021–2100 Random Forest layer .................................... 69
Figure 35. PHM in Madii Lii for SSP585:2021–2100 Random Forest layer ................................ 70
Figure 36. British Columbia testing and training points for BIOCLIM SDM. ............................. 71
Figure 37. British Columbia BIOCLIM layer (SSP126:2021–2040). ........................................... 72
Figure 38. British Columbia PHM MaxEnt layer for SSP126:2021–2040. .................................. 74
Figure 39. British Columbia PHM Random Forest layer for SSP126:2021–2040. ...................... 75
Figure 40. British Columbia PHM Random Forest layer for SSP245:2021–2100. ...................... 77
Figure 41. British Columbia PHM Random Forest layer SSP585:2021–2100. ............................ 78
Figure 42. British Columbia BIOCLIM for (SSP126:2021–2040). .............................................. 90

vi

Acknowledgement / Dedication
I would like to express my appreciation to Luutkudziiwus community for allowing me onto
their traditional territory and to undertake this research. In particular, I would like to thank Wingchief Richard Wright who sadly passed away in May of 2020. Richard’s mission was to reclaim
control of his, and his community's traditional territory and to gather data on the severity and
longevity of damage exacted on Lax’yip Madii Lii by the Crown and its agents over the last one
hundred and fifty years. As Pansy Wright-Simms has now succeeded Richard as Wing-Chief,
she has echoed the voice of her brother and endorsed the continuation of this research. Thank
you.
I would also like to sincerely thank Dr. Chelsey Geralda Armstrong for her supervision and
support during this project. Her patience and guidance have been monumental. She has been an
inspiration and allowed me to become a better student and person. With Dr. Armstrong’s support I
was awarded the Natural Sciences and Engineering Research Council of Canada, Canada Graduate
Scholarships—Master’s program award, for which I am so incredibly grateful.
I would also like to thank Kenny Rabnett, Madii Lii natural resources manager, who provided
direction and the background for the project, and Nabby who provided incredible food and
accommodation for the field data collection. I also thank Allie Golt for being a brilliant research
assistant during the field data collection.
I would like to recognize the efforts of Dr. Roger Wheate, who replaced my previous cosupervisor whose abrupt departure left me somewhat orphaned, alongside Dr. Phil Burton who joined
my committee. Dr Burton’s knowledge, of and experience with devil’s club has been essential for the
completion of this project. I sincerely acknowledge the assistance and direction provided by Dr. Alex
McAlvay for his guidance through the conception stage and for providing support with the species
distribution modelling section of this research. I would also thank my fellow graduate students,
Andrew Boxwell, for being there in times of need, and all the team at UNBC’s Applied Analysis
Hub, Sunny Tseng, Lisa Koetke and Juan Delgado. Finally, I thank my amazing wife, Rachel, for her
patience, understanding and help along this winding road.

vii

Chapter 1
1.0 INTRODUCTION
The intergenerational transfer of cultural and spiritual knowledge about traditional medicinal
plants is at grave risk due to deforestation and the increasing threat of climate change
(Applequist et al., 2020; Buenz, 2005; Sharma et al., 2020). Millions of Indigenous peoples
worldwide use traditional herbal medicines as their principal form of healthcare, and as a way of
supporting their daily livelihoods. Up to 80% of the world’s population depend on medicinal
plants for healthcare and income generation (Farnsworth & Soejarto, 1991; Uprety et al., 2012;
Zhang et al,, 2002).
Since the imposition of colonial states and attendant resource extraction practices,
Indigenous communities have witnessed dramatic changes in their landscapes. Globally,
Indigenous peoples own and occupy over a third of the earth’s remaining intact forested
landscape (Fa et al., 2020). In regions where Indigenous peoples continue to maintain and
steward their lands, forest biodiversity and productivity are overwhelmingly higher than in areas
managed by state governments (Dawson et al., 2021; Garnett et al., 2018). For example, areas of
rainforest cultivated and managed by Kayapo people in Brazil’s Amazon Basin are more
productive and have superior biodiversity value compared with the surrounding forest (Posey,
1984). Similarly, throughout the world, the use of fire stewardship in forested ecosystems by
Indigenous peoples has maintained habitat heterogeneity and promoted biodiversity, whereas
forested land that has undergone colonial/western-style fire suppression is now experiencing
frequent abnormally severe wildfires, negatively impacting biodiversity (Bowman & Haberle,
2010; Hoffman et al., 2021; McWethy et al., 2019).Unfortunately, the practice of unrestricted
legal and illegal deforestation continues to increase in order to support the expansion of

1

agriculture, mining, and oil and gas extraction. For example, in the Brazilian Amazon, annual
deforestation dramatically decreased by 80% between 2005 and 2012, but since that low-point,
deforestation in the area has increased by 122% per year in 2019 compared with 2012 (Carrero
et al., 2020), significantly contributing to greenhouse gas emissions (GHGs) and disrupting
global precipitation (Gedney & Valdes, 2000). Regrettably, similar stories of deforestation
directly impacting Indigenous peoples are replicated worldwide and continue to threaten
communities with alienation from their traditional lands. In many of these cases, Indigenous
communities are at the forefront of maintaining impactful stewardship strategies, central to both
long- and short-term climate change mitigation (Masson-Delmotte et al., 2021).
The global demand for timber and timber products has increased steadily, with a
significant acceleration after World War II (McNeill & Engelke, 2016; Warman, 2014). In
Canada, the timber harvesting industry underwent a period of rapid modernization in the 1950s.
New mechanized equipment, such as skidders (high-clearance four-wheel drive tractors) and
powerful diesel-fueled log-transport trucks, made previously inaccessible forest-stands possible
to log. Coupled with the new equipment, the highly efficient method of harvesting known as
clearcutting became popular. “Clearcutting” is described as “the removal of the entire stand in
one cutting, with reproduction obtained artificially (by planting or seeding) or by natural
seeding from adjacent stands or trees cut in the clearing operation” (Keenan & Kimmins, 1993,
p. 122). In British Columbia (BC), clearcutting became common-place for the logging industry
during the 1960s, and has left an indelible mark on the province that remains evident today.
Industrial tree-harvesting in northwestern BC has caused a myriad of detrimental impacts
to forested ecosystems, such as slope failure and mass wasting events. For example, in Rennell
Sound, Haida Gwaii, research showed that where logging and accompanying forest road

2

construction had occurred, there was a 34 x increase in landslide frequency and 31 x increase in
sediment yield on disturbed areas (Schwab, 1983). Similarly, in 1994, 2.5 million m3 of debris
was released during the Mink Creek landslide, near the town of Terrace. The failure occurred at
the site of two clearcut harvesting events, one in 1981 and the second in 1985. The failure was
likely caused by a combination of successive years of elevated precipitation and the liquefiable
nature of the glaciomarine silty-clay sediments. However, the research also showed that almost
half of the landslide area had been clearcut and that the retrogression was arrested at the interface
with the intact mature forest (Geertsema et al., 2006).
In the Lakelse Lake area of the Skeena Region, the construction of over 300 km of
forest access roads between 1964 and 1972 caused a dramatic increase in bedload sediments
entering the river system, with associated channel instability, and damage to salmonid habitat
(Gottesfeld et al, 2002).] Studies showed that 43 of the 64 stream reaches associated with the
Lakelse watershed were highly impacted by logging-related sedimentation and had suffered a
reduction in productive capacity (Rabnett, 2006).
For decades, government decision-makers and policymakers in British Columbia have
argued that highly regulated forestry practices are sustainable and that they replicate the impacts
of natural disturbances, or that silvicultural practices can lead to similar or even accelerate stand
development (Kramkowski, 2012). In one example, Barker et al. (2014) argued that disturbances
like fire and infestation have similar effects to clearcut logging, using the instance of the McLure
Fire of August 2003 in central British Columbia. The study compared seedling establishment and
growth performance in five disturbed areas, classified as High Severity Burn, Low Severity
Burn, Clearcut, Screefed Clearcut, and Undisturbed Forest. The study found seedling growth in
the Screefed Clearcut was comparable to the High Severity Burn treatment. Therefore, logging

3

can reproduce comparable effects to intensive wildfire regarding seedling establishment, growth
and nutrition (Barker, et al. 2014). Increasingly, however, scientists recognize that these studies
could be more robust and acknowledge the important role of forest dynamics such as soil
nutrient change acidification and compaction and that critical temporal constraints are often
ignored (Haeussler & Kneeshaw, 2003). When comparing clearcuts with intensive fires, standlevel disturbances are complex and can differ widely depending on biogeoclimatic zone, variant ,
elevation and aspect, and can be altered by additional factors such as pine beetle infestation
(McIntire, 2001) and canopy-gap size (Coates & Burton, 1997). For the Interior Cedar-Hemlock
biogeoclimatic zone, age-class data shows that fire return ranges from 101 years in the dry
ICHdw subzone to 129 years in the ICHmw2 variant, to 226 years in the wet ICHwk1 variant
(Pollack et al. 1997). However in the Arrow Forest District in southeast BC, disturbance
estimated intervals vary dramatically within both the ICHdw subzone and the ICHmw2 variant,
with ranges of 65–449 years and 97–458 years, respectively (Wong et al., 2003).
Industrial-scale logging in the Skeena Region of BC has generated economic growth for
the region, and eventually led to the construction of the Columbia Cellulose pulp mill, which
precipitated the expansion of the logging industry in the region (Wedley, 1990). From the 1950s
onwards, the wage-labour economy played a progressively significant role for Indigenous
communities and sustained local economies (Menzies & Butler, 2008). Many members of the
Ts’msyen and Gitxsan communities participated in forestry in the early industrial period. Some
worked for Columbia Cellulose, which held the first Tree Farm Licence (TFL) and logged
Nisga'a and Kitsumkalum territories during the 1960s and 1970s (Menzies & Butler, 2008).
Nevertheless, debates continue to emerge about the industry’s long-term viability, particularly

4

given the disparate wealth gaps it has created, and the boom-bust cycling of local economies and
environmental impacts that are increasingly felt in an age of climate change.
To regulate industrial-scale logging, Forest Stewardship Plans (FSPs) are required prior
to timber harvesting, and are meant to assess logging impacts on a range of values, including
biodiversity, fish and riparian areas, soils, visual quality, water quality and recreation (BC
Ministry of Forests and Range, 2009). Similarly, heritage impact assessments, legislated under
the Heritage Conservation Act (1996), assess logging impacts on cultural and archaeological
heritage, particularly Indigenous peoples’ cultural heritage. Unfortunately, these pieces of
legislation fail to address concerns and issues that arise from logging in areas where there are
culturally valued plant species (Lepofsky et al., 2020). In Gitxsan contexts, concerns over the
impacts of logging on devil's club (Oplopanax horridus) populations are not considered under
any legislation or regulation. This is a troublesome gap, given the extreme value and salience of
this plant for Gitxsan people, and the recognition that local access to desirable patches is
diminishing.
In order to bridge this gap, this research attempts to quantify and qualify the potential
impacts of logging on devil's club populations. Specifically, the purpose is to (1) bridge
regulatory gaps in the impact assessment process, which divides cultural-ecological values,
therefore omitting important components of logging impacts on Gitxsan health, well-being,
cultural practices, and livelihoods; and (2) take a historical-ecological approach to assess logging
impacts on devil's club over time. Temporal scales, as a component of cumulative effects, are
often downplayed in the context of impact studies. As such, I use a chrono-sequence approach
centred within the Luutkudziiwus (Gitxsan) homeland of Madii Lii.

5

1.1 Background
1.1.1. Ethnobotanical Context
Devil’s club, known as wa'uumst or hu'ums in Gitxsanimx or Gitksenimx, is a perennial,
shade-tolerant, understory deciduous shrub that grows throughout northwestern North America,
from lowland to subalpine zones, within coastal and mountainous temperate coniferous forests. It
is generally found in greater profusion in regions with higher precipitation and on benches,
alluvial fans, stream edges and seepages, and flourishes in calcareous soils (Alaback, 1984;
Klinkenberg, 2021). Devil’s club typically occupies canopy gaps in temperate forests with welldrained, moist, and nitrogen-rich soils (Hemstrom & Logan, 1986; Klinkenberg, 2021). This
species occupies various environments and elevations from sea level to 2750 m, and often
favours steep, north-facing slopes (Klinkenberg, 2021).
Devil’s club produces long spiny stems, which can grow to 3 m tall but can grow up to 5
m in length, often developing along the forest floor in a recumbent position (Lantz, 2001; Lantz
& Antos 2002). The stems are protected with thick, yellow spines 5–10 mm long, which can
cause infection if left untreated (Drummond, 2018; Klinkenberg, 2021). The plant produces
palmately lobed leaves, with between seven to nine lobes per leaf, ranging in size from 10 cm to
35 cm wide with a heart-shaped base (Klinkenberg, 2021). In the late spring the racemes
elongate, up to 25 cm in length, and are covered with short-stalked, greenish-white flowers
(Figure 1). Later in the summer, the flowers are replaced with bright red berries measuring 5–8
mm in length (Figure 2) (Klinkenberg, 2021).

6

Figure 1. Devil’s club green/white flower, budding in early June. Photo by Haeussler S. (2023)

Figure 2. Devil's club red berries ripening in mid-August. Photo by Smith A. (2021)
As with many understory species with limited resources, vegetative reproduction is the
primary source of expansion, and in this case devil’s club, it uses clonal layering to reproduce.
7

Within a patch, new shoots or ramets attempt to occupy gaps to access sunlight, creating a
network of crossing stems. Over time the combined weight of the stems forces the bottom-most
stems into the soil and expands the patch (Pickett & Kempf, 1980; Roorbach, 1999). Eventually,
these ramets produce their own root system and detach from the host plant, becoming a separate
entity. In addition to vegetative reproduction, devil’s club also possesses the ability to reproduce
sexually, though it appears to be infrequent, and it is challenging to quantify how often it is
successful (Kachel et al., 2006; Lantz 2001; Lantz & Antos 2002; Roorbach 1999).
Devil’s club is used and named by 38 linguistic groups across northwestern North
America (Turner, 2014). The healing properties of devil’s club have likely been employed by
Indigenous peoples over millennia, and remain one of the most commonly used medicinal plant
species (Calway et al., 2012; Compton, 1993; Turner, 1982, 2014).
An early account of the medicinal use of devil’s club was documented by the chief
physician of the Russian-American Company, Eduardo Blaschke in 1842, who observed Tlingit
using it to treat sores (Blaschke, 1842 fide Lantz et al., 2004). Decades later, anthropologist
Harlan Smith, employed by the Canada Department of Mines and the National Museum of
Canada, recorded Indigenous communities using devil’s club for various medicinal purposes. His
documentary films showed how members of Bella Coola (Nuxalk) Nation chewed devil’s club
bark to treat whooping cough (Smith, 1927b). Dakelh (previously Carrier) drank a tea made from
crushed inner-bark as a purgative and as a treatment for dysentery (Smith, 1927a). Smith also
recorded how Gitxsan Elders made a decoction of boiled bark and twigs as a remedy for
consumption (tuberculosis) (Smith, 1927c).
In 1982, renowned ethnobotanist Nancy Turner and her Indigenous colleagues catalogued
the multiple ways devil’s club is used by Indigenous peoples of the Pacific Northwest (Turner,

8

1982). Lantz et al., (2004) showed the variety of preparations derived from devil’s club as
treatments or medicines for more than 34 types of sicknesses or conditions, as well as many
spiritual applications. Many of the treatments use a combination of inner bark, root, and root
bark to treat arthritis, rheumatism, and diabetes. It is also used as a purgative, laxative, and
analgesic for pneumonia, coughs and colds and used for sores, swellings, cuts, bruises and skin
infections (Lantz et al., 2004; Turner, 2014). For Gitxsan community members, devil’s club is
used to treat tuberculosis, gynecological issues, and venereal disease (Gottesfeld & Anderson,
1988a; Johnson, 2006a, 2013; Turner, 2014). [there should be semi colon between Johnson,
2006a and 2013 – they are separate references]. The use of devil’s club as a treatment for
diabetes is commonly cited by Gitxsan and neighbouring Ts’msyen Elders (Johnson, 2006). The
Northern Health Authority in British Columbia has begun experimenting with trials to support
Elders in taking devil’s club for diabetes. Over recent years there has been a burgeoning interest
in the healing properties of devil’s club and salves manufactured from devil’s club under the
commercially adopted names of “Alaskan Ginseng” and “Pacific Ginseng” (Lantz et al., 2004;
Tilford, 1997; Vuksan et al., 2019).
Multinational pharmaceutical corporations have explored the phytochemicals extracted
from devil’s club to develop potential treatments for conditions such as cancer and diabetes
(Calway et al., 2012; Cheung et al., 2015, 2019; Tai et al., 2014; Wang et al., 2013). Some of the
most significant advances have been made in the treatment of colonic carcinoma cell apoptosis
(Wang et al., 2019), pancreatic and ovarian cancers, as well as a potential treatment for Type II
diabetes (Cheung et al., 2019; Elahmer, 2018; Yang et al., 2018).
In addition to the medicinal benefits and values of devil’s club, the act of harvesting is
recognized as a spiritually sanctioned practiced. Conventionally, only those considered to be

9

specialist medicinal practitioners are authorized to harvest devil’s club, and generally prefer to
harvest alone and in remote locations (Turner & Berkes 2006; Turner et al., 2004). Protocols for
devil's club harvesting are complex and differ for different communities, however the concepts
of sustainability and reciprocity are repeated throughout the Ts’msyen and Gitxsan canon. For
example, Gitga’at (Ts’msyen) protocols note that no more than four stalks should be harvested at
a time and the harvester must offer a prayer or leave a gift of tobacco or a coin in the spot from
where the plant was taken (Turner & Berkes 2006; Turner et al., 2004). In Haida culture, the
numbers four and forty reoccur often and are regarded as having “power”. For example, Haida
fishermen and hunters would perform a complex ritual for good luck, taking 40 sticks of devil’s
club; each piece measured from the elbow to the fingertips. These sticks were pushed into the
sand to form a circle, and the fisherman/hunter would sit inside. Next, they would use their teeth
to gnaw four girdling strips of the inner bark from each of the 40 sticks and slowly chew and
swallow them. The ritual would induce nausea and vomiting, eventually leading to a form of
hallucination or shamanistic vision, bringing good luck to the trip, with the additional benefit of
making the subject attractive to women (Curtis, 1916 fide Deagle, 1988).
Previous research indicates that the preferred (size) diameter of devil's club stems for
medical or spiritual reasons is larger than 1" diameter (25 mm) (Burton & Burton, 2015). In
addition to the significance of stem diameter, the harvest location is a central consideration
(Johnson, 2006; Johnson, 2013). It is said that devil’s club from higher (montane) elevations has
a greater potency than that harvested at a lower elevation, because the brief growing season
concentrates the active chemical compounds and the plant benefits from more intense ultra-violet
light at higher altitudes (Larson, 1988; Turner et al., 2011).

10

1.1.2. Logging Context
Throughout British Columbia, forestry has dramatically altered the land cover, and land-use
practices since the surge of settler colonial activity, beginning in the mid 19th century. For
example, prior to 1827, the forests in the Lower Fraser Basin consisted primarily of coniferous
stands on uplands and lowlands in the eastern portion of the basin (North & Teversham, 1984)
By 1930, urbanization and agriculture had taken over much of the delta area and lowland near
the Fraser River, and approximately 50% of the region’s forests were cleared for agriculture.
Most of the Sitka spruce (Picea sitchensis) and yellow cedar (Callitropsis nootkatensis) were
removed (McCombs & Chittenden, 1990), while Douglas-fir (Pseudotsuga menziesii), western
redcedar (Thuja plicata), and western hemlock (Tsuga heterophylla) comprised most of the
remaining timber . Douglas-fir was common at low elevations, and western redcedar and
hemlock was common at higher elevations (British Columbia Department of Lands, 1930). Over
time, the less accessible cedar-hemlock forests became logged. By 1990, the amount of
coniferous forest in the Fraser Basin had increased due to the reforestation policies implemented
during the 1930s. However, the overall forested area (including deciduous forest) was
significantly lower than before settler colonialism. Any remaining mature forest was primarily
located at higher elevations or areas with challenging access (Boyle et al., 1997; McCombs &
Chittenden, 1990). This sequence, of logging floodplains, valley bottoms, benches and lowerelevation shallow slopes, was replicated throughout the province, including in the Skeena Region
of northwestern BC. Furthermore, by the 1960s, advancements in tree harvesting technology
made accessing highly profitable cedar-hemlock timber from steeper and higher elevation terrain
possible. Combined with the new practice of clearcut harvesting increasing threats to devil’s club
throughout British Columbia has been a concern.

11

Recently, forestry experts and resource managers have begun to recognize that the
impacts of forestry may threaten devil’s club habitat. For example, Canfor’s Sustainable Forest
Management Plan (SFMP) for the Prince George Region, states that ecosystems such as “Devil’s
club - Ostrich fern” are currently threatened ecosystems and should be protected during future
forest harvesting (Kaye, 2012). In addition, the provincial government recognizes there are 8
blue-listed (considered to be of special concern) and 2 red-listed (candidates for extirpated,
endangered, or threatened status) devil’s club based ecological communities that are either rare
or at risk (BC Conservation Data Centre, 2023). For example, red-list status has been designated
for Thuja plicata / Oplopanax horridus / Matteuccia struthiopteris found in the ICHvk2/05
Biogeoclimatic Unit, which has been incorporated into the Provincial Forest and Range Practices
Act (FRPA). The Act is intended to minimize the effects of forest and range practices on
threatened species or habitat and specify a commitment to maintaining their current and historic
ranges. Likewise, Picea engelmannii x glauca / Oplopanax horridus / Hylocomium splendens
ecosystem found in the Biogeoclimatic Unit, SBSmc1/07 is also red-listed, but is yet to receive
FRPA protection.
A study undertaken in the Lakelse-Kitimat River Valley over from 1986 to 1987
examined the effects of both natural disturbances and forestry on devil’s club and thimbleberry
(Rubus parviflorus). The study found that both thimbleberry and devil’s club grew relatively
better when canopy cover decreased due to disturbance. The study found that thimbleberry
apportioned more of its resources to flowering within canopy-gaps, whereas devil’s club showed
an increase in leaf size. The results suggested that, as an early-succession species, thimbleberry
was able to out-compete devil’s club and better establish in the high-light disturbed habitat,

12

whereas devil’s club was more slow-growing and arrived considerably later in the succession
stage (Mason, 1990).
Other research in central British Columbia indicated that twenty years after clearcut
logging and post-harvest slash burning, many understory plant species were able to recover.
However, several key species, including white-flowered rhododendron (Rhododendron
albiflorum) and devil’s club, had yet to re-establish over the period of study (Chandler et al.,
2017). Similarly, research undertaken in the Sitka spruce-western hemlock forests in southern
Alaska found that salmonberry (Rubus spectabilis), Alaskan blueberry (Vaccinium ovalifolium)
and red huckleberry (Vaccinium parvifolium) recovered within 30 years post-clearcutting.
However, devil’s club was only present in significant quantities in the mature even-aged forest
between 100 and 250 years old (Alaback, 1984).
To date, there are only a handful of studies documenting the extent of timber harvest
impacts on ethnobotanically important species for Gitxsan people, including devil’s club.
Previous research about the effects of clearcut logging on devil’s club in British Columbia
indicates a complex interplay of factors that affect the persistence and regrowth of devil's club
after timber harvesting. During their research on Vancouver Island, Lantz & Antos, (2002) found
that devil's club can persist through clonal expansion following clearcutting. However, there was
significantly increased ramet mortality, and as could be expected, growth was slow in logged
areas compared to unlogged areas. In a more recent study also located in northern Vancouver
Island, Beese and colleagues (2022), observed that in four different silvicultural systems
(clearcutting, patch cutting, dispersed retention, and shelterwood) adjacent to old-growth forest,
the only shrub to exhibit continued adverse effects from logging was devil’s club. For their pilot
study in the Nass Valley near Terrace, Burton and Burton (2015) noted that the thickest stemmed
13

devil’s club (over 25 mm diameter), which are preferred for traditional medicinal use, were
uncommon and difficult to find in replanted forest-stands that were less than eleven years old.
Further research in the Date Creek Experimental Forest in northwestern British Columbia
indicated that devil’s club had been negatively impacted by clearcut logging, during both initial
logging events and during forest recovery (Lilles et al., 2018).
These findings leave outstanding questions about the impacts of logging on devil’s
club—particularly how populations are impacted and over what temporal range. Furthermore,
most studies assessing the impact of logging on devil’s club is relatively short-term (2–20 years).
Lastly, BC is relatively culturally and biologically diverse, impacts studied on Vancouver Island,
for example, do not necessarily scale up to studies in Gitxsan territories.
1.2 Research Objectives
This study aims to understand the effects of timber harvest on devil’s club in Lax’yip Madii Lii,
one of two territories of Wilp (house) Luutkudziiwus (Gitxsan), between Smithers and Hazelton,
along the Suskwa River drainage in northwestern British Columbia. Concerns over logging in the
region are not new. Objections to illegal logging and land theft in Gitxsan territories led to the
precedent setting Delgamuukw and Gisday'wa v. the Crown (S.C.C.44) in the Supreme Court of
Canada. For millennia, Gitxsan house groups, including Wilp Luutkudziiwus, have successfully
managed their valued resources by increasing productivity of desired habitats for animals and
plants, for food and medicinal purposes (Armstrong, 2017). However, clearcutting accelerated in
the 1960s and continues to have compounding effects and impacts on the Madii Lii landscape.
The near extirpation of the high-value western redcedar from Madii Lii, and the ongoing
fragmentation of the Interior Cedar-Hemlock (ICH) biogeoclimatic zone (a preferred habitat for
devil's club), has drastically reduced the principal environment where devil's club can flourish.
14

Today, the landscape is a patchwork of former cutblocks with immature and young pine
forest dominating stands logged in 1970s and 1980s, through to grassy scrub communities, with
sparsely scattered shrubs and herb assemblages from those sites logged in the late 1990s and
2000s. The logging industry may have reduced availability and possibly disrupted access (even
with the construction of logging roads) to preferred devil's club harvesting sites, therefore
impeding the transfer of significant cultural knowledge and practices concerning this valuable
medicinal plant.
To assess logging impacts on devil’s club in Madii Lii, I analyze devil’s club
distributions in relation to timber harvesting events and consider the health and vigour of patches
throughout in the territory. Overall, the goal is to consider the role of logging history in shaping
the structure of devil’s club habitat, compared with sites that have been maintained by
Luutkudziiwus, and without industrial logging.
With their determination to assert Aboriginal rights and title, promote a cultural
resurgence, and engage with younger community members, in the last decade Luutkudziiwus
leaders and house members constructed a culture camp, also called Madii Lii. The camp fosters
pride in Gitxsan history and promotes values, language, and cultural knowledge, and provides
opportunities for members to re-engage in traditional practices and livelihoods such as fishing,
hunting, and medicinal and food plant harvesting—including the harvest and preparation of
devil’s club. Without access to high-quality (e.g., thick-diameter stems, isolated stands) this
fundamental knowledge is consistently challenged. Therefore, plant morphology data collected
during the research process could be used as an inventory of the locations of the longest and
thickest stemmed devil's club specimens in Madii Lii. In addition, these data make it possible to
calculate the time interval required after logging before devil's club reaches the desired size for
15

medicinal and cultural use to be harvested; therefore, a harvesting schedule and land-use plan
that consider important plants like devil’s club could be generated to ensure future availability.
This research was designed to help inform Wilp Luutkudziiwus land-use planning
objectives by identifying important aspects of devil’s club habitat, distribution, and access. For
example, data collected on stem length and diameter — both important characteristics for harvest
— were assessed, mapped, and disseminated to house members. I also use species distribution
models to assess where populations currently thrive, should thrive, and if/how those aspects of
habitat determination will change under various climate change scenarios.
Table 1. Thesis research objectives and the research questions that need to be answered to
achieve those objectives.
Research Objective
Research Question
a) What is the impact on devil’s club caused by
1) Compile a series of recommendations for
logging on the culturally salient medicinal plant
future forest management practices and
for wilp Luutkudziiwus?
landscape-level planning objectives
b) Where are the locations of the most accessible
concerning the health and location of
devil’s club plants of a preferred size suitable for
culturally preferred devil’s club specimens
cultural purposes?

2) Understand the local variation in the

c) What environmental factors, such as slope,

presence/absence of devil’s club within

aspect, and elevation, play a role in the size

Madii Lii

of devil’s club?
d) What is the nature of the current access and
availability of devil’s club in the territory?

3) Understand the effects of logging on
the abundance and morphology of
devil’s club in Madii Lii.

e) Compared to unlogged areas, is there a
relationship between devil’s club stem
abundance, length and diameter in forest

16

harvest cutblocks and the time since logging
(TSL)?
4) Examine the effect of shade on the

f) Is there a relationship between leaf area and

abundance and morphology of devil’s
club in Madii Lii.

shade provided by the forest canopy?
g) Is there a relationship between shade and
devil’s club leaf area, on devil’s club stem
abundance, length and diameter in forest
harvest cutblocks from different TSL
compared to the unlogged areas?

5) Use species distribution models to

h) What is the current species distribution for

understand the projected impact of
future climate change on devil’s club’s
presence

within

Madii

Lii

devil’s club in Madii Lii?
i) What is the expected species distribution for

and

throughout British Columbia.

devil’s club based on climate change
prediction

models

in

Madii

Lii

throughout British Columbia?
1.3 Species Distribution Models Literature Review
Species distribution models, also known as environmental (or ecological) niche models
(ENMs), or predictive habitat models (PHMs), are used in biology, conservation ecology and
quantitative ecology to project expected variations in the geographic ranges of species as the
effects of global warming become increasingly apparent (Iverson & Prasad, 2001; Lawler et al.,
2009; Thuiller, 2004; Thuiller et al., 2005). (Peterson et al., 2011; Peterson & Soberón, 2012).
For example, SDMs have been useful for predicting potential extirpations and/or extinctions of
species over broad time scales (Condro et al., 2021; Franklin, 2010; Gregory et al., 2012). By
combining environmental, climate, and geospatial data describing the current range and habitat
associations of a species, an ecological envelope is defined, and the probability of it persisting

17

and

within a distinct multidimensional environmental envelope can be determined (Bunce et al.,
2013; Guisan & Zimmermann, 2000; Wong et al., 2007).
Conservation biologists and resource managers often use SDMs to determine the range of
a species of particular interest, for example plant species that are used for food, technology, or
medicinal purposes. In Egypt, a study used the SDM algorithm MaxEnt (Phillips et al., 2006) to
generate species-richness maps for 114 medicinal plants growing in protected areas (Kaky &
Gilbert, 2016). The SDM data showed that species richness was significantly higher in protected
areas than outside, and conservation management—at least for those species—was effective
(Kaky & Gilbert, 2016). In the Indian state of Jharkhand, researchers used SDMs to understand
the impact of over-harvesting and changes in land-use and land cover on the deciduous tree
species Buchanania cochinchinensis (Mishra et al., 2021). Indigenous peoples of the region use
virtually every part of this plant for various medicinal purposes, including treatment of fever,
ulcers, diarrhea, dysentery, asthma and snakebite (Siddiqui et al., 2014). The research determined
where the region’s most suitable habitat is today and identified potentially suitable sites on which
to focus conservation planning in the future (Mishra et al., 2021). The pharmacologically and
commercially valuable plant species American ginseng (Panax quinquefolius) has been grown
commercially in China since the 1980s. However, cropping has become increasingly difficult in
recent years due to continuous over-harvesting, soil degradation, and reduced genetic diversity
(Cruse‐Sanders and Hamrick, 2004). Therefore, researchers have generated SDMs using 37
ecological variables based on precipitation–temperature ratios to establish a suitable model to
determine the location of potential new producing areas of American ginseng using the
maximum entropy model. (Zhang et al., 2018).

18

A handful of studies have utilized SDMs to understand the impact of deforestation on
culturally important and medicinal plant species (Benner et al., 2019; Nimasow et al., 2016; Wan
et al., 2016). A study based in British Columbia's Great Bear Rainforest used SDMs to help
predict the distribution of “monumental’ western redcedar (Benner et al., 2019). The research
showed that this culturally important tree was disproportionally affected by industrial logging
and negatively impacted the traditional livelihoods of Indigenous communities in the region
(Benner et al., 2019). Monumental redcedar are typically over 250 years old, and are
exceptionally large, high-quality trees suitable for cultural purposes such as carving dug-out
canoes and totem poles, and for building long houses (e.g., Garibaldi & Turner, 2004).
Recently the most common application for SDMs is for predicting future distributions of
taxa which are confronted with the increasingly dire consequences of climate change.
Researchers have used SDMs to determine the expected range of medical and culturally
important plants growing within the Pacific Northwest (Bogawski et al., 2019; Hirabayashi et al.,
2022; Mucioki et al., 2022; Prevéy et al., 2020; Prevéy et al. 2020b). For example, Prevéy et al.,
(2020a; 2020b) examined the potential impacts of climate change on three culturally relevant
plant species (used for food and medicine) for Indigenous peoples across the Pacific Northwest,
including beaked hazelnut (Corylus cornuta), Oregon grape (Mahonia aquifolium), and salal
(Gaultheria shallon) In another study, researchers examined how change in weather patterns,
precipitation, and temperature will affect the phenology of huckleberry (Vaccinium
membranaceum) and how this will impact the range of this vital food source (Prevéy et al.
2020b).
Finally, with the simplification of SDMs and the perception that SDMs are now
generating more accurate species distribution predictions, they have become increasingly
19

popular. Such modelling, provides conservation biologists with a valuable tool for planning the
management of endangered species (Olden et al., 2010). However, some scholars have indicated
that there are fundamental concerns about data availability and quality (Dormann, 2007; Santini
et al., 2021). This suggests there could be a danger in over-relying on SDM when making critical
decisions and shaping policy (Sinclair et al., 2010). All SDM research should be done critically
and with all limitations honestly and succinctly laid out.

20

Chapter 2
2.0 RESEARCH METHODS
2.1 Introduction
To understand the effects of industrialized logging on devil’s club over five decades, I used
geospatial and statistical data analysis, and predictive habitat modelling methods based on fielddata collected in Madii Lii in 2021 (Figure 5). I compared devil’s club growing in timber
harvested areas from decadal periods, beginning in the mid-1970s onwards until the last largescale logging events in 2012. Timber harvest sites with devil’s club were then compared with
devil’s club patches growing in relatively undisturbed sites (controls) and sites that were utilized
by Gitxsan people, but had not been logged (e.g., cultural use sites, archaeological sites, etc.).
In order to achieve these objectives, I employed a four-phase approach. First, I
characterized potential devil’s club habitat to ascertain potential sampling locations. Second, I
conducted fieldwork in the territory, searching potential sampling locations. Where devil’s club
was present, I collected habitat data and individual plant morphology measurements. Third,
habitat and morphology data were analyzed to determine any correlations or patterns using the R
Studio data analysis platform. Finally, devil's club species distribution models were constructed
to map potential distribution under current habitat conditions and future climate scenarios.

21

2.2 Study Area
Madii Lii lies approximately 300 km inland from the British Columbian coastal port-city of
Prince Rupert, within the Babine sub-basin watershed of the Skeena River, northeast of the town
of Hazelton (Figure 3).

Figure 3. Madii Lii and the Northwest coast of British Columbia.
The Suskwa River watershed defines most of the boarders around Madii Lii territory, which
extends as far as the Suskwa-Bulkley Rivers confluence to the west, and the Babine drainage to
the northeast. The highest elevation landform in the territory is found Hahla Gyoot (Thoen
Mountain) at 2285 m in the Babine Range. The lowest point of elevation in the territory is where
Fifteen-Mile Creek joins the Suskwa River at 360 m, with the mean basin elevation being 1160
m (Rabnett, 2018). The climate in Madii Lii is largely transitional, being warmer and wetter than
the interior climate to the east but cooler and drier than the coastal climate to the west
(Meidinger & Pojar, 1991). The territory is characterized mostly by the Interior Cedar-Hemlock
22

(ICH) biogeoclimatic zone (variants moist-cold-1 and moist-cold-2), representing the transitional
influence of maritime and continental climates at lower elevations. Engelmann Spruce-Subalpine
Fir Biogeoclimatic zone (ESSF) is strictly interior (continental) and reflects the mountainous
terrain and shorter growing seasons found at higher elevations. The upper limits of the territory
are characterized by treeless alpine tundra of the Boreal Altai Fescue Alpine (BAFA) Zone,
which has more boreal (northern interior) affinities than coastal/southern interior mountains
(Banner et al. 1993).

Figure 4. Luutkudziiwus’ Xsi Gwin Hauums and Madii Lii territories.
Study sites were sampled within two ICH variants present in the territory: the Nass Moist Cold
(mc1) variant and the Hazelton Moist Cold (mc2) variant. The dominant tree species in the ICH

23

zone are western hemlock (Tsuga heterophylla), spruce (Picea), lodgepole pine (Pinus contorta
var. latifolia), and subalpine fir (Abies lasiocarpa). The ICHmc1 variant lies slightly higher than
the ICHmc2, between 750 m and 1100 m. In the ICHmc1, seral species such as trembling aspen
(Populus tremuloide) and paper birch (Betula papyrifera) are often found on southern aspects.
Step moss (Hylocomium splendens), red-stemmed feathermoss (Pleurozium schreben), and
knight's plume (Ptilium crista-castrensis) cover the forest floor, and the herb and shrub-layers
comprise of black huckleberry, oval-leaved blueberry (Vaccinium ovalifolium), and false azalea
(Menziesia ferruginea) (Banner et al. , 1993; Burton et al., 2000). The ICHmc2 subzone
appears at 350 m lower elevation and consistently ascends to 750 m. At the lower elevations
human activities occur more frequently. It differs mainly from the mc1 subzone due to the
presence of small, scattered stands of western redcedar and also contains a much greater
abundance of broad-leaved deciduous trees such as aspen, paper birch and black cottonwood.
Similar to the mc1 variant, the most abundant ground cover is step moss, with herbs and shrub
layers consisting of huckleberry and blueberry (Banner et al., 1993). The ICHmc1 variant blends
into the Engelmann Spruce-Subalpine Fir (ESSF) zone at higher elevations (900 m-1300 m
above sea level), depending on localized variations in topography, aspect, and microclimates.
Subalpine fir and western hemlock co-dominate in this zone with lodgepole pine (Pinus
contorta), and Engelmann spruce (Picea engelmannii) interspersed throughout. The ground
cover comprises of a mixture of red-stemmed feathermoss and leafy liverworts (Barbilophozia
spp.). The shrub layer consists of black huckleberry and false azalea, with a variety of ferns, and
five-leave bramble (Rubus pedatus) is a prominent species in the herb layer (Banner et al., 1993,
Burton et al., 2000). The Boreal Altai Fescue Alpine (BAFA) zone is found at higher elevations
above the timberline and supports small amounts of ground vegetation. Madii Lii’s climate and

24

weather patterns are typically produced in the northeastern Pacific Ocean from the “Aleutian
Low” This extratropical-low typically produces severe winter storms, often containing
exceptionally high-moisture content, which precipitates as either rain or snow depending on local
temperatures. The Suskwa Valley weather station (Environment Canada-107G879) is located
adjacent to the Suskwa-Harold Price confluence at an elevation of 534 m. The 35-year record
shows mean annual precipitation is 757 mm, of which 551 mm is rainfall and 206 mm as snow
(Environment Canada 2010). Temperature measurements indicate the average minimum winterlow as -9.3oC, average winter-low of -6.4oC and the maximum summer-high of 20.8oC with the
average summer-high of 14.5oC over the last 35 years (Environment Canada-107G879).

2.3 Site Selection
Madii Lii territory covers a 352 km2 area, encompassing mountainous terrain and deeply incised
valleys. There are over seven decades of logging history in Madii Lii territory, which has
resulted in a patchwork mosaic of regenerating cutblocks at different stages of succession and
forest recovery. To subdivide the territory into manageable sampling sites representative of the
various stages and temporal scales, and to increase the likelihood of locating devil’s club
patches, I compiled remote sensing data, habitat data, previous ecological studies in the region,
and timber harvest data. A table was created to characterize accessible devil’s club habitat based
on features such as slope, elevation, aspect, and logging history (Table 2). The table was
compiled using data from United States Geological Survey (USGS) and mapping data acquired
from the Province of British Columbia Ministry of Forest, Land, Natural Resources Operations
and control sampling sites.

25

Table 2. Ecological criteria, parameters and analysis method to determine the location of control
sampling sites.
Variable

Parameters

Analysis Method

Data source

References

Not Clearcut Logged

Sites that have not
experienced clear cut
logged

Visual assessment of
satellite image

LandSat_8 (2018-0727) WRS PATH 52
ROW 21

USGS satellite image
Landsat 8

Forest type

Coniferous forest

Visual assessment of
satellite image

LandSat_8 (2018-0727) WRS PATH 52
ROW 21

USGS satellite image
Landsat 8

Forest type

Coniferous forest

vector analysis

British Columbia
Vegetation Resource
Inventory

(British Columbia.
FLNROD 2018a)

Age of stand

Lead species over 100
years old

vector analysis

British Columbia
Consolidated Harvest
Areas

British Columbia.
FLNROD 2018b)

Dominant ICH tree
species

Lead species is typically
found in ICH mc1/mc2

vector analysis

British Columbia
Forest Tenure License

(British Columbia.
FLNROD 2017)

ICHmc1 elevation

750 - 1100m

raster analysis

Digital Elevation
Model (DEM) T093M25m resolution

(Canada 2019)

ICHmc2 elevation

350m-750m

raster analysis

Digital Elevation
Model (DEM) T093M25m resolution

(Canada 2019)

Slope

slope between 110 and
350 (avg 210)

raster analysis

Digital Elevation
Model (DEM) T093M25m resolution

(Canada 2019)

North aspect

3550 - 00

raster analysis

Digital Elevation
Model (DEM) T093M25m resolution

(Canada 2019)

Northeast aspect

00 – 450

raster analysis

Digital Elevation
Model (DEM) T093M25m resolution

(Canada 2019)

Water bodies

Within 50 m of moving
or static water bod

vector analysis

British Columbia
Rivers streams, lakes
wetlands atlas

British Columbia.
FLNROD 2018d)
(Roorbach 1999)

located within 1000 m
of Arterial FSR

Forest Service Roads
(FSR)

vector analysis

British Columbia
Forest roads atlas

(British Columbia.
FLNROD 2018c)

located within 500 m
of active/deactivated
resource roads

Open or abandoned
forest roads

vector analysis

British Columbia
Forest roads atlas

(British Columbia.
FLNROD 2018c)

located within 500m of
trail system

Gitxsan Trails System

vector analysis

Gitxsan Cultural
Heritage Trail System

(Rabnett 2018)

Edge effect mitigation

120 m wide buffer zone
around internal
perimeter of potential
control site

vector analysis

120m no-sample edge
zone within potential
devil's club control
sites

(Harper et al., 2005)

26

To further characterize and predict the location of typical devil’s club habitat, supporting
information was obtained from a variety of sources, including Biogeoclimatic Ecosystem
Classification zones (BEC) program (BC Ministry of Forests and Range 2021). In addition, the
vegetation tables from Banner et al. (1993) describe where and on which site-series devil’s club
is likely to occur in each biogeoclimatic variant and in what relative abundance. I also assessed
research such as MacKinnon et al. (1992), Lantz and Antos (2002), Burton and Burton (2015),
and Mackenzie and Meidinger (2018) to further refine my search areas for representative devil’s
club habitat. These locations meeting the criteria outlined in Table 2 were then overlayed with
the Kispiox Predictive Habitat Mapping (PHM) used for grizzly bear habitat modelling, which
provided a high-resolution, site-series description of the habitat and ecology of the region
(Hanchard, 2021). These GIS data were then processed to identify potential devil’s club habitat
and sampling sites within the ICH mc1 and ICH mc2 BEC zone using a weighted overlay
(ArcMap, 2021). Because access to field sites containing devil’s club was essential for the
research, but also for Luutkudziiwus community members for their continued access to high
quality harvesting opportunities, the distance from Gitxsan trail system (500 m buffer) and forest
roads (1000 m buffer) was used as a matrix parameter to determine potential field-sampling
locations.
In order to assess the effects of logging on devil's club habitat, timber harvest data from
Harvested Areas of BC (FLNROD 2018b) was used to build a geodatabase of timber harvest
events in the territory, what we refer to as Time Since Logging (TSL). To compare areas that
have been logged with areas that have not (i.e., controls) we referenced the Vegetation Resources
Inventory (VRI) Province of British Columbia (FLNROD 2018a). We selected each sample site
based on logging events over four ten-year period intervals: 1970–1979, 1980–1989, 1990–1999,

27

and 2000–2012. Throughout the subsequent chapters, I refer to these as “groups”, that is, Group
1 (1970–1979), Group 2 (1980–1989), Group 3 (1990–1999), and Group 4 (2000–2012), (Figure
5).

Figure 5. Devil's club sampling sites (24) on Luutkudziiwus traditional territory Madii Lii.

All preferred habitat data and timber harvest history data were then overlayed to
recommend potential sampling sites and control sites. It was apparent that many of the potential
control sites were adjacent to timber harvested areas, so a buffer zone of 120 m was established
around all potential control sites to reduce the influence of any edge effects (e.g., higher light or
wind due to nearby canopy removal, drainage impacts, etc.).

28

Fifty potential sample sites were identified using ArcMap 10.8.1 as either control sites or
logged sites, including 8 control sites, and 42 logged sites from 1970–1979 (n=12), 1980–1989
(n= 12), 1990–1999 (n=12), and 2000–2012 (n=6). In addition, I selected four (n=4) samples
sites that were neither logged, nor technically “controls” but had undergone land-use changes by
Wilp members over the last 200+ years. These “traditional use” sites are located near to the
Suskwa River, along the ancient Gitxsan trail system, now known as The Babine Trail, as well as
near archaeological village sites and sites with important place names, such as Nox, Nox Lobbit
(“Faces of Our Ancestors”).

2.4 Data Collection
Field data were collected in June 2021. A visual assessment of the pre-selected sites (50 sites)
was undertaken to confirm the accessibility to the site and the presence/absence of devil’s club.
In total, 24 suitable sites (26 sites were unreachable due to access issues) were sampled (Figure
5), including: 4 control sites, 4 traditional-use sites, and 16 timber harvest sites from 1970–1979
(n=5), 1980–1989 (n= 4), 1990–1999 (n=4), and 2000–2012 (n=3). The logged stands sampled
ranged between approximately 20 ha and 60 ha with a mean area of about 40 ha, with forest
stands identified as control sites and traditional use sites being considerably smaller. At each site,
I recorded the location (at the mid-point of the site perimeter, at the 30 m mark perpendicular to
the forest road) and elevation, and surveyed the area for devil’s club (presence/absence). Tree
cover at each site was described by percent cover of tree species in the overstory and was
estimated visually. A general description of the site was noted (evidence of burning, proximity to
water, understory shrubs), but was not used in analyses. At each site, four circular plots were
sampled at 20 m intervals on an 80 m transect (Appendix 1). Each plot had a 5.64 m radius,

29

equating to a total sampling area of 100 m2 plot (0.01 ha). Where devil’s club was absent, a
maximum of two plots were sampled. At each plot I recorded (1) a range of ecological variables,
(2) devil's club presence/absence and cover, and (3) in each plot a series of morphometric
measurements were taken as a proxy for devil’s club health and vigour.
2.4.1. Ecological Variables
At each plot, I characterized ecological and geophysical characteristics, including 1) forest type,
(e.g., coniferous, deciduous or open), (2) dominant tree species in descending order of frequency,
(3) slope, (4) aspect, (5) canopy closure, and (6) understory plant assemblages. The slope was
obtained using the Measure® software application on an Apple iPhone. Each plot’s location,
elevation; and aspect were recorded by using the applicable functions available on a Garmin
handheld GPS. Canopy openness was measured using the %Cover™ (version 2.1) canopy
application on an Apple iPhone. The phone/camera was placed at the centre of the plot at
shoulder height, and 4 photos were taken from different directions and a mean value of percent
canopy cover was determined. The software immediately processed the image for light and shade
and then expressed the amount as a percentage, and it was then recorded onto the datasheet. To
describe the relative ground cover for understory species of shrubs, herbs, mosses and lichens, I
used the Cover Abundance Scale (Blanquet 1964) method. The Cover Abundance Scale uses a
1–5 ranking system, where : 5 = >75 percent cover; 4 = 50-75 percent cover; 3 = 25-50 percent
cover; 2 = 5-25 percent cover; 1 = numerous, but less than 5 percent cover, or scattered, with
cover up to 5 percent. Any unidentified specimens were photographed for later identification,
and any presence of water, either flowing or standing water were recorded.

30

2.4.2. Presence-Absence Data
At each plot, devil’s club presence-absence and stem abundance was recorded. Only devil’s club
stems greater than 30 cm in length (including recumbent stems) were recorded, because the stem
diameter measurement was taken 30 cm above the stem base. Therefore, any specimen plant <
30 cm height would not be able to have its diameter recorded. However, to ensure that all
presence-absence data was recorded, any stems < 30 cm were noted separately and recorded and
included in presence-absence modelling
2.4.3. Devil’s Club Measurements
In each plot, a series of morphometric measurements were taken as a proxy for devil’s club
health and vigour (Dawe et al., 2017). New growth was measured from the leaf scars defining
the end of the previous year’s growth to the tip of the terminal bud using a flexible centimetre
measuring tape. The stem diameter was measured at 30 cm above the base of the stem with
digital-calipers. A representative sample of four leaves was removed from each plot, then
labelled and measured for leaf area. The four representative leaf samples were obtained by taking
a single large leaf from a plant in each quadrant of the plot (wherever possible), and if there was
not an available plant in each quadrant, the next-nearest plant was selected. All the leaf samples
which had been collected and labelled earlier in the day were processed at the fieldwork
accommodation camp in the evening. The samples are positioned flat on a whiteboard, and leaf
area was then calculated using the iPhone application LeafByte™ (version 1.3.0) (GetmanPickering et al., 2020). Four dots surround the leaf in a square; these constitute the scale for area
calculations; a 30 cm square frame was used for most of the larger leaves and 20 cm for the
smaller samples. The four dots gave a frame of reference for determining absolute areas in the
image and for image-size correction if the image was taken at an angle. A photograph was taken,

31

and the image was processed to provide leaf area in cm2, and the results were recorded in a
spreadsheet for later analysis.
2.5 Data Analysis
Site-level data (slope, elevation, tree cover, plant assemblages) were analyzed and assessed along
with plot data (devil's club stem length, stem diameter, stem abundance, leaf area, and shade),
and site use history (e.g., control, traditional use, or timber harvest site). The statistical analysis
software environment R Studio Version 4.1.1 (R Core Team, 2021) was used for all data
analysis. A Shapiro-Wilks (1965) test was employed to examine the distribution normality of
stem abundance, stem length and stem diameter.
Due to the non-normal distribution of the data, a non-parametric significance test
(Kruskal-Wallis, 1952) was used to evaluate the effect of time since logging (TSL) or standage/time since disturbance (TSD) on plot-level data. Then a posthoc Dunn (1964) test was used
to test for any significant relationships between individual pairs (e.g., control sites with Group 1).
I then used the Bonferroni-Holm method (Holms, 1979) p-value adjustment to adjust for multiple
comparisons. To illustrate the effect of site history on stem abundance, stem length, and stem
diameter, a series of boxplots were generated. Next, a Spearman Rank Correlation Test package
(Savicky, 2022) was used to identify any relationship between pairs of these three variables.
It was expected that site-elevation would have a confounding effect on the data with
regards to devil’s club presence/absence, stem abundance, stem length, and stem diameter when
comparing sites which experienced different histories –Time Since Logging (TSL) and stand
age/Time Since Disturbance (TSD). Therefore, testing was performed on the TSL (industrialscale logging) and stand age/TSD for control and traditional use sites while controlling for the
32

influence of elevation on the morphometric variables. First, the normality of distribution was
tested using the Shapiro-Wilks test. To examine the relationships between TSL and elevation on
stem abundance, stem length, and stem diameter, a Kendall-Theil Sen Siegel non-parametric
analysis of covariance (Siegel, 1982) was employed in conjunction with the R Studio Median
Based Linear Models (mblm) package (Komsta, 2013). The mblm analysis provides linear
models based on Theil-Sen single median, and Siegel repeated medians. This formula provides a
robust 29 or 50 percent breakdown point, respectively. For each point, the slopes between it and
the others are calculated (resulting n-1 slopes) and the median is taken. This results in (n)
medians and median from this median then calculates the slope estimator. Consequently, the
only additional values provided are slope and the intercept, and it does not provide a Multiple Rsquared value. Therefore, to attain an R2 value a pseudo-R-squared value was generated with the
R-Studio script “efronRSquared” (Mangiafico & Mangiafico, 2017)
Finally, given devil’s club proclivity for shaded environments, I tested stem abundance,
stem length, stem diameter, and size of individual leaves as relates to tree cover by using
Spearman correlation. Due to the non-normal distribution of the data, any regression analysis
required non-parametric tests. Kendall-Theil Sen Siegel non-parametric analysis of covariance
was performed for each variable. Finally, the frequency of devil’s club occurrences was
analyzed with respect to the presence of the differing tree species associated with each plot and
was illustrated with a bar chart.
2.6 Species Distribution Modelling
The final component of this research was to determine current devil’s club habitat availability in
Madii Lii and for its potential distribution across the province of British Columbia, and to assess

33

how habitats may be altered due to impacts from climate change. To do this, I used the presenceabsence data from the fieldwork to create a species distribution model for Madii Lii and data
from the online repository Global Biodiversity Information Facility (GBIF) (www.gbif.org,
2022) to create a species distribution models for devil’s club throughout British Columbia. In
addition, all modelling bioclimatic data was retrieved from the WorldClim® (worldclim.org)
global data-repository (Fick & Hijmans, 2017). To generate the SDMs, a workflow employing
computational operations was performed in the R version 4.1.1 environment (R Core Team
2021) in a stepwise procedure that incorporated several statistical analysis packages. In
preparation for the modelling, a map representing Madii Lii and devil’s club presence-absence
data-points recorded during the data collection was generated in ArcMap 10.8.1 (ESRI, 2020)
(Figure 6.). The SDM modelling package “dismo” version 1.3.5 (Hijmans & Elith, 2021) was
used as the foundation for model development on the R Studio platform. The “dismo” (version
1.3.5) package applied several species distribution model profiles, BIOCLIM (Nix, 1986) and
Domain (Hijmans, 1996) models, and also provided a link from the R environment to the

Figure 6.Devil’s club species occurrence within Madii Lii recorded during fieldwork data
collection in June 2021.

34

maximum entropy ‘MaxEnt’ (version 3.4.3) model (Phillips et al 2006) General Linear
Regression/Gaussian (GLM) model (McCullagh & Nelder, 1989), and the Support Vector
Machine (SVM), (Vapnik, 1998) models. Hence, all five SDM model profiles were used during
the model generation process. The stepwise process was as follows: (1) BIOCLIM analysis used
the 17 bioclimatic variables derived from WorldClim 2 (www.worldclim.org) model for the
current period (1970-2000) to derive the entire dataset to generate the initial base SDM, and was
visualized with seventeen individual graphic outputs displaying each variable
(www.worldclim.org; Fick & Hijmans, 2017) (Appendix 2). (2) Having generated the multiple
BIOCLIM outputs, the data were then used to extract raster values and generate data points
assigned to the bioclimatic variables as predictors. (3) The model algorithm then used the
BIOCLIM layers to generate testing and training points from the presence-absence data collected
during fieldwork, and these points were then used to evaluate the model-quality by producing a
Relative Operating Characteristics (ROC) plot.
The Relative Operating Characteristics (ROC) plot indicates the quality of the model,
initially for the BIOCLIM, then for subsequent model profiles. The model output was assessed
using the Area Under Curve (AUC) for the receiver operator curve (ROC). The AUC is a
measure of rank correlation. The higher AUC indicates that sites with high predicted suitability
values are areas of known presence. A satisfactory AUC value would range from 0.6 to 0.7 for
the model. When AUC is approximately 0.5, the model has no discrimination capacity to
distinguish between positive class and negative class. The higher the AUC, the better the model
performs at differentiating between the positive and negative classes. (4) A response plot was
produced to show the relationship between the probability of the occurrence of devil’s club with
35

each environmental variable. The plots then indicated which variable had the most influence on
devil’s club distribution. (5) The Domain, GLM, MaxEnt and SVM models were tested, and their
AUC output values were compared with each other to ascertain which was the best performing
model (Appendix 4). (6) The “average score model” was produced by combining the output from
BIOCLIM, Domain, GLM, and SVM. In addition, a weighted mean of three models (BIOCLIM,
GLM and SVM) was also generated but was not used (Appendix 5).
2.6.1. Madii Lii Predicted Habitat Model
In order to estimate the future availability of devil’s club in Madii Lii for Luutkudziiwus
community members, I generated species distribution models for two different climate scenarios.
This analysis used the Shared Socioeconomic Pathways (SSPs) scenarios (Riahi et al., 2017;
IPPC AR6 2021) to illustrate how global warming and the subsequent disruption to localized
weather patterns may affect devil's club viability within the territory and the species range
provincially. The analysis used the Future climate projection from Canadian Earth System Model
version 5 (CanESM5) Swart et al., (2019) in conjunction with the Shared Socioeconomic
Pathway (SSP) as the basis for climate modeling. The analysis uses five different scenarios,
SSP1–SSP5, based on projected global socioeconomic changes (used as a predictor for
variability in greenhouse gas emissions) from the present day through to the end of the 21st
century. SSP126 represents current conditions in 2022 with approximately 420 ppm atmospheric
carbon, which was used as the baseline for the current conditions and an increase of 1.7o C by
2040. For the predictive habitat modelling, I used the scenarios SSP245, intermediate GHG
emissions: CO2 emissions remain around current levels until 2050, then falling but not reaching
net zero by 2100, with an expected increase in temperature ranging from 2.0o C to 3.5o C by the
end of the century. For the worst-case scenario I selected SSP585, very high GHG emissions:

36

CO2 emissions triple by 2075, which leads to predicted global temperature increasing from 3.3o
C to 5.7o C (CarbonBiref.org 2018). By using the species distribution models generated from the
previous section and then employing future scenario simulations SSP245 and SSP585, I provided
predicted future habitat suitability for devil’s club. The presence-absence data collected during
the fieldwork in June 2021 were used to construct the models for Madii Lii and presence data
obtained from GBIF repository for British Columbia (GBIF.org, 2022) with the current condition
SDM, the BIOCLIM model was used to generate the first model and allowed me to evaluate
which of the 19 variables contributed the most influence over the presence of devil’s club. The
model was then assessed using the AUC value. Using the data obtained from the BIOCLIM
model, I then ran the Random Forest (Breiman, 2001) and MaxEnt models, and performed
evaluation using the AUC score to determine which model performed best. Random Forest was
then selected to produce the predictive SDM modelling for the climate scenarios. I then created
predictive habitat suitability models for both (SSP245 and SSP585) for each 20-year period,
2021–2040, 2041–2060, 2061–2080 and 2081–2100 using the Random Forest model.

Figure 7. British Columbia recorded locations of devil’s club within the province (www.gbif.org, 2022).

37

2.6.2. British Columbia Predictive Habitat Model
Species distribution and habitat suitability models were generated to illustrate where devil’s club
can be expected to persist and/or change distributions within the Province of British Columbia
under two predicted climate change scenarios. As with the climate change predictive models
generated for Madii Lii, I employed the moderate carbon increase level scenario SSP245 and the
worst-case scenario model SSP585. First, I defined the British Columbia boundary in ArcMap
10.8.1 and obtained devil’s club presence data from the Global Biodiversity Information Facility
(GBIF) (www.gbif.org) (submitted by iNaturalist membership, and concentrated around
population centers) within that geographical territory (3642 occurrences within North America).
I then clipped the species occurrence data to include only presence data that occurred within
British Columbia (1600 occurrences) (Figure 7). The 19 bioclimatic variables were employed to
generate the testing and training model and to generate the dataset using the k-fold crossvalidation. For the devil’s club predictive model to be accepted as being applicable, a rigorous kfold test was performed. This technique ensured that each observation from the training and
testing dataset which determined the presence or pseudo-absence (otherwise known as
background data points) had a reasonable probability of appearing in the k-fold cross-validation.
For this scenario, 5-fold cross-validation (k=5) was used. The test and train dataset were then
applied to the base SSP126 (2021–40) climate scenario using the BIOCLIM profile to produce a
base model for current conditions for British Columbia. The next operation was to generate the
Maxent environmental bioclimatic variable contribution table/graph, and the derived
environmental bioclimatic variable response plots for the Maxent predictive habitat model for
SSP126 (2021–40). The MaxEnt model was then evaluated for accuracy using the AUC output
value.

38

To ensure that the research produced the most accurate predictive model possible, both
the Random Forest (Breiman, 2001), and Support Vector Machine (Vapnik, 1998) profiles were
tested and compared using their AUC scores as indictors for best model fit. As Random Forest
returned the best performing model, it was then used to complete the predictive habitat suitability
models for both scenarios. Predictive habitat suitability models were generated for both (SSP245
and SSP585) for each 20-year period, 2021–2040, 2041–2060, 2061–2080 and 2081–2100 using
the Random Forest model profile.

39

Chapter 3
3.0 RESULTS
In this chapter, results are summarized in two sections. In the first, I present results of field data
collection and analyses, including comparisons of devil’s club habitat and morphology across
sites with different land-use histories. Second, I present a series of predictive habitat models for
devil’s club in Madii Lii and across British Columbia. These models are meant to evaluate where
suitable devil’s club habitat will persist and occur under various climate change scenarios at both
a regional and continental scale.

3.1 Devil’s Club Morphometric Analysis
Devil’s club was present in all three site types: traditional use sites, timber harvest sites, and
unlogged control sites. The Shapiro-Wilks test reveal that values for devil’s stem abundance,
stem length and stem diameter at all sites were not normally distributed (Appendix 3). As nonparametric analyses were required for these data, the one-way Kruskal-Wallis test was combined
with the Dunn post-hoc test for multiple comparisons of groups. Furthermore, the Kendall-Theil
Sen Siegel non-parametric linear regression models were favoured for evaluating devil’s club
stem abundance, stem length and stem diameter throughout the analysis.

3.1.1. Abundance, Stem Length, Stem Diameter
The results shows there is a significant difference (p = 0.0019 H-test= 18.985, df = 5) in devil’s
club stem abundance across the three sampled site types (timber harvest, traditional use, and
controls). Devil’s club stem abundance was not significantly different across the five decades of
timber harvest sites. The plot shows there is fewer stems in the 1970s sites than in the 1980s then
there was a gradual decrease in the abundance of stems from timber harvested sites in most

40

recently logged cutblocks in the 2000s. The Dunn's post-hoc test indicates devil's club stem
abundance is significantly greater (Table 3) in traditional use sites than in all timber harvest sites.
There is no significant difference between the control sites and timber harvest sites, but there is
between the control sites and traditional use sites (Table 3). The traditional use sites group
showed a median abundance of 15 stems/100 m2 plot, which is considerably more than that of
the harvested and control areas (Figure 8).

Figure 8. Distribution of stem abundance per 100 m2 (0.01 ha) devil’s club at each type of site.
The bold line within the boxplot represents the median value of the variable. The box represents
the interquartile range. The lines extending parallel from the interquartile range. The lines
extending parallel from the boxes are known as the “whiskers,” which are used to indicate 1.5 x
the 25th or 75th percentile values.

41

Table 3. Pairwise results from Dunn post-hoc tests (with Bonferroni adjusted p-values)
comparing devil’s club stem density at traditional use sites with harvested sites and control sites.
Time periods in parentheses indicate the period of logging.
Comparison

Z - statistic

P. unadjusted

P. adjusted

Control - Traditional use

-3.0101

0.0026

0.0339

G1 - Traditional use (1970–1979)

-3.3179

0.0009

0.0127

G2 - Traditional use (1980–1989)

-2.0449

0.0409

0.1170

G3 - Traditional use (1990–1999)

-2.3607

0.0182

0.2030

G4 - Traditional use (2000–2012)

-3.6950

0.0002

0.0032

.

The Kruskal-Wallis test reveals that devil’s club stem length is significantly different (p=0.0155
H = 14.019, df = 5) among all sampled sites (timber harvest, traditional use, controls).
When comparing between the unlogged and logged sites stem length differs most significantly
between the most recently logged sites, 2000–2012 and the traditional use sites (p=0.0282), and
when comparing between logged-sites from different decades the most significant difference is
between Group 2, logged 1980–1989 and Group 4, logged 2000–2012 (p = 0.0061), where in both
cases the former was significantly longer than the latter. Generally, the median devil's club stem
length from the earlier logged sites is greater than from the more recently logged sites, except for
Group 2 (1980–1989) stems, which were longer that Group 1 (1970–1979) stems (Figure 9).

42

Figure 9. Distribution of median stem length of devil’s club plants in different stand ages of
forest harvest, control group and traditional-use sites.

Data analysis results from the Kruskal-Wallis test reveal that devil’s club stem diameter is
significantly different (p= 0.0102, H= 15.036, df = 5) among all sampled sites (timber harvest,
traditional use, controls). Though not significant, the analysis shows that stem diameter differs
the most between traditional use sites and the control sites (p=0.0504), closely followed by the
traditional use sites and Group 4 (2000–2012) (p=0.0566). The largest median stem diameter is
found on the traditional use sites (16.4 mm), and the smallest median diameter from those sites
harvested from 2000–2012 (12.4 mm) (Figure 10).

43

Figure 10. Distribution of median stem diameter of devil’s club plants found growing in different
stand ages after forest harvest, including control sites and traditional use sites.
3.1.2. Examine Relationship Between Shade and Leaf Area on Stem Length
The results of the Kendall-Theil Sen Siegel non-parametric linear regression indicates that both
shade (p<0.0001) (residual standard error: 46.78 on df= 59), and leaf area (p< 0.0001) (residual
standard error: 117.3 on df=59) shows a significant positive correlation between stem length and
shade (Figure 11). The scatter plot indicates that as the shade increases, the stem length also
increases. Also, a positive correlation can be seen between stem length and leaf area. This shows
that as the leaf area increases, the stem length also increases (Figure 12).

44

Figure 11. Relationship between devil’s club stem length and canopy closure.

Figure 12. Relationship between devil’s club stem length and devil’s club leaf area. Dashed line
indicates a relationship but not a cause and effect.
.
45

3.1.3. Examine Relationship Between Shade and Leaf Area on Stem Diameter
The results from the Kendall-Theil Sen Siegel non-parametric linear regression show that there is
significant relationship between shade and devil’s club stem diameter (p = 0.0391; residual
standard error: 3.866 on df=59). The analysis indicates an overall negative correlation between
shade and stem diameter; as shade increases, stem diameter decreases. The scatter plot, shows
that many of the smallest diameter stems are growing in >75% shade (Figure 13). Conversely,
there is a significant positive association of leaf area and stem diameter (p<0.0001; residual
standard error: 10.34 on 59 degrees of freedom). The analysis indicates a relationship between
leaf area and stem diameter, such that when the leaf-area increases, the stem diameter also
increases (Figure 14).

Figure 13. Regression analysis indicating relationship between devil’s club stem diameter and
canopy shade.

46

Figure 14. Relationship between devil’s club stem diameter and devil’s club leaf area. Dashed
line indicates a relationship but not a cause and effect.
.
3.1.4. Relationships Between Devil’s Club Metrics
Across all sites, the Spearman's rank correlation test reveals a strong correlation between stem
diameter and stem-length, indicating that as devil’s club stems increase in length, they have a
relatively proportionate increase in their basal diameter (p-value < 0.0001, ρ=0.6436). In
addition, the relationship between stem diameter and stem abundance also shows a significant
correlation (p-value< 0.0001 ρ = 0.4678), and likewise stem length and stem abundance show a
significant positive correlation (p-value < 0.0001 ρ= 0.4396).

47

3.1.5. Relationship Between Elevation and Devil’s Club Metrics
In order to control for the variation in elevation across all sampled sites, a Kendall-Theil Sen
Siegel regression analysis was used to examine devil’s club stem abundance, stem length, and
stem diameter as related to site elevation. The regression analysis shows that stem abundance
had a significant relationship with elevation across all sites (p<0.0049; residual standard error:
7.318 on df=59). The analysis indicates that as the elevation increased, the stem abundance
decreased, and consequently, the most productive areas were found at lower elevations,
approximately 400 m to 500 m above sea level, which corresponded with the location of the
traditional use sites situated in the Suskwa River drainage (Figure 15)

Figure 15. Relationship of devil’s club stem abundance (stems/100 m2 plot) and elevation using
Kendall-Theil Sen Siegel regression analysis.
The regression analysis show that stem length has no significant relationship with elevation
across all sites (p= 0.698; residual standard error: 46.96 on df=59). Though a linear trend was
not significant, many of the sites with the longest stems are found at mid-elevations,
48

approximately 600 m to 800 m, above sea level (Figure 16), which corresponded to harvested
sites from the earliest harvested periods during the 1970s and also the control sites.

Figure 16. No significant relationship of devil’s club stem length (cm) and elevation using
Kendall-Theil Sen Siegel regression analysis (p>0.05).
Regression analysis further indicates that stem diameter has a significant relationship with
elevation across all sites (p=0.002, residual standard error: 3.575 on df=59). The analysis
indicates that as the elevation increases, the stem diameter decreases, and overall, the thickest
stems were found at lower elevations, many of which were located in the traditional-use area,
with numerous patches containing more than 20 stems over 25 mm diameter. In contrast, at
higher elevations, above 900m, there are few sparsely scattered, large diameter stems. However,
several large, thick-stemmed specimens found in sites located between 600 m and 800 m
elevations corresponding with sites logged during the 1970s and 1980s (Figure 17).

49

Figure 17. Relationship of devil's club stem diameter (mm) and elevation (m using Kendall-Theil Sen
Siegel regression analysis).
When comparing devil's club stem abundance, stem length and stem diameter with site history,
Time Since Logging for stands < 60 years old or Stand Age/Time Since Disturbance for those
areas traditionally logged or initiated by natural disturbance (e.g., fire) for all sites while
controlling for the influence of elevation, the results show a significant relationship between
stem abundance and TSL (p = 0.0009; residual standard error: 7.805 on 59 degrees of freedom;
Figure 19). The results for stem length show no significant relationship between the TSD across
all sites, even when the effect of elevation was removed (p= 0.1700; residual standard error:
46.19 on 59 degrees of freedom) (Figure 18). However, for stem diameter with TSL there is a
significant relationship between (p=0.0002; residual standard error: 4.015 on 59 degrees of
freedom). The results from the analysis show that as TSL increased, both stem abundance
(Figure 19) and stem diameter increase (Figure 20).

50

Figure 18. Regression analysis plot shows the relationship between devil’s club stem length and
the Stand Age or Time Since Disturbance (TSD)Time Since Logging while controlling for the
influence of elevation (p>0.05).

Figure 19. Regression analysis plot shows the relationship between devil's club stem abundance
and the Stand Age or Time Since Disturbance (TSD)using a multi-regression to control for the
influence of elevation.

51

Figure 20. Relationship of devil's club stem diameter and the Stand Age or Time Since
Disturbance (TSD) using a multi-regression to control for the influence of elevation.
3.2 Canopy Openness and Leaf Area Analysis
The removal of tree canopy from timber harvesting may or may not affect the size of devil’s club
leaves (leaf area), which in turn may affect their health and vigour (Dawe et al., 2017). To test
the relationship between canopy openness and leaf area across all three stand history categories, I
use a Spearman’s rank correlation test. Results indicate that, across all sites there is no
significant relationship between canopy openness and the area of individual devil’s club leaves
(p = 0.1412, ρ=0.1845). However, evaluating leaf area across all sites a Kruskal-Wallis nonparametric test show that there was a significant relationship between leaf area. And TSL
(p=0.0028 H = 18.111, df = 5). Specifically, these data reveal that median leaf size was greatest
on traditional use sites, with the oldest timber harvest site, Group 1 (1970-1979) having larger
leaf area than the rest of the timber harvest sites. These data also indicate that the specimens with

52

the largest median leaf area are in areas which has the most time to regenerate since logging
events (Figure 21).

Figure 21. Distribution of the median leaf area (cm2) of devil’s club plants found growing in
different stand ages of forest harvest.

3.2.1. Effect of Shade and Leaf Area on Abundance
The results of the Kendall-Theil Sen Siegel non-parametric linear regression model demonstrates
that shade does not have a significant relationship on devil’s club stem abundance (stems/100 m2
plot (p=0.0566, residual standard error: 8.299 on 86 degrees of freedom; Figure 22), on the othe
hand there is a significant correlation from leaf-area on the abundance of devil’s club stems
(stem/100 m2 plot) (p<0.0001, residual standard error: 20.78 on df= 86). The analysis shows a
positive relationship between leaf area and abundance. The regression line (Figure 23) suggests
that on sites where the leaf area increases, devil’s club stem abundance also increases.

53

Figure 22. Relationship between devil’s club stem abundance and available shade (p>0.05) using
regression analysis.

Figure 23. Relationship between devil’s club stem abundance and devil’s club leaf area. Dashed
line indicates a relationship but not a cause and effect.

54

.
3.3 Tree Cover analysis
The relationship of devil’s club abundance to tree cover is presented in the bar chart below (Figure
24). Of the 10 tree species encountered mountain alder had the lowest number of occurrences at
locations where devil’s club was present (n=6). In contrast, the most frequently occurring tree
species on sites where devil’s club was present (n=50) was spruce (various natural hybrids) . The
analysis showed that across all sites, spruce (Picea spp.) (30.1%) was the most frequently
occurring species where devil’s club was present, followed by western redcedar (15.7 %) and paper
birch (10.2 %).

Figure 24. Frequency of occurrences of tree species where devil’s club was present.

55

3.4 Species Distribution Modelling Results
The development of the species distribution model for devil’s club in Madii Lii successfully
followed the stepwise workflow as described in the methods section of this research. The
BIOCLIM profile used the 17 bioclimatic variables to produce the base model and determine
which variable had the most influence over the occurrence of devil’s club. These data were then
also used to generate the devil’s club testing and training data points (Figure 25). For model
evaluation calculations the BIOCLIM algorithm generated 18 species presence data-points and
144 species absence data-points from the 88 presence-absence data-points (plots) obtained
during the field data collection (Table 4). The model also generated “background points” which
were created to characterize the environmental qualities of the study region. The model
algorithm was then able to determine which conditions the target species (devil’s club) would
most likely need to persist without manufacturing absence locations.

Figure 25. Training and testing points generated to produce a BIOCLIM species distribution
model. The blue points represent devil’s club present training locations and the green points
denote devil’s club absence locations, yellow are background testing.

56

3.4.1. Response Plot
The results from the bioclimatic variables are illustrated as response plots for each individual
variable. The response plot for devil’s club provides the range of values in each predictor
variable based on the fitted models for the SDM. By using the presence-absence data gathered
during the fieldwork, a bioclimatic response graph can be generated for each variable. Each plot
indicates the relationship between the probability of occurrence for devil’s club with respect to
each individual environmental variable (Figure 26).

57

Figure 26. Multiple response plots for devil’s club probability of occurrence in relationship to all bioclimatic
variables evaluated on the basis of training data. The x-axis represents the range of values of the environmental
variable, and the y-axis provides the probability of incidence on a scale from 0 (low probability) to 1 (high
probability) (Worldclim, 2021).
58

3.4.2. Interpretation of Response Plots
The results from the generated response plots showed the relationship between the probability of
occurrence of devil’s club relating to each environmental variable. For each response plot, the
output was generated for one environmental variable. The response plots for Bio(s)
1,5,6,8,9,10,11,14,15 all displayed similar shaped profiles. The response curves for these plots
showed that the probability of occurrence of this species follows an optimum curve. For example
the response curve for Bio 1 Mean Annual Temperature, indicates a high probability of species
presence in areas that maintain a temperature between ~0 and 4 oC. The response plot curves for
the Bio(s) 4, 12,16, and 17 indicate the probability of species presence over a wide range of
values. For example, Bio 12, Annual Precipitation, the curve shows that devil’s club prefers
wetter climates up to a maximum of 730 mm per annum, then as precipitation becomes too
excessive the probability of occurrence declines. Finally, the response plot curves indicate that
for Bio(s) 2, 3, and 7, all have a very steep and narrow response mean, indicating that devil’s
club has a small tolerance for variations in these variables, for example, Mean Diurnal Range,
the mean of monthly (max temp–min temp).
3.4.3. BIOCLIM Model Evaluation
This BIOCLIM model indicated an AUC value of 0.8951 which suggested a very reliable model.
The AUC of 0. 8951 value translates to 89.51% chance that the model will reliably differentiate
between positive and negative classes (Table 4). In addition, the model evaluation relies on two
variables, (p) presence and (a) absence, and uses these values to calculate the correlation
coefficient value (0.453). Furthermore, the max TPR+TNR value also provides an alternate
model evaluation. The True Positive Rate (sensitivity) and the True Negative Rate (specificity)
across the range of the threshold probability values can be calculated/plotted to provide a
complementary model evaluation. The value represents the probability threshold at which our

59

model maximizes the True Positive Rate (TPR) and the True Negative Rate (TNR). The greater
the value for both figures, the likelihood of a more reliable model
Table 4. Summary of BIOCLIM model evaluation
Class
Model Evaluation
n presences
18
n absences
144
AUC
0.895
correlation coefficient
0.453
max TPR+TNR
0.293

The interpretation of the predicted map is quantified as the suitability of a grid cell (n-row 24, ncolumn 68, n-cell 1632, 6889 m2 cell) on a scale from 0 to 1, where 0 refers to very low
suitability, and one refers to very high suitability. The model based on BIOCLIM layers
combined to generate a predicted presence-absence for devil’s club in Madii Lii (Figure 27).
Figure 27. Left: species distribution model with continuous habitat suitability values determined
from the BIOCLIM layer model. Right: binary presence/absence model used by applying a
threshold, derived from the max TPR+TNR value.

3.4.4. General Linear Regression Method for the SDM
The General Linear Regression/Gaussian model results produced another set of predicted habitat
suitability and presence-absence maps (Figure 28). The suitability map shows current conditions,
with the green (darker) area where the most suitable habitat for devil’s club is more likely to be

60

found. The yellow areas are moderately suitable, the pink-orange are low suitability, and the
beige/white are non-suitable habitats for devil’s club. The presence-absence diagram highlights
in green the region to the south where devil’s club is expected to proliferate.

Figure 28. Left: species distribution model with continuous habitat suitability values for the Generalized Liner
Model/Gaussian model. Right: binary presence/absence model used by applying a threshold, derived from the max
TPR+TNR value.
3.4.5. Average Score Model
The process of model generation was then completed for the Domain, MaxEnt and SVM models
and the AUC scores compared for likelihood of best fit for each model (Appendix 4). By
merging the results of four models (BIOCLIM, Domain, GLM and SVM; Figure 29) an “average
model” was calculated by employing a mean score of predicted probability of devil’s club
occurrence from each of the individual contributors to provide a single average model prediction
(Figure 30).

61

Figure 29. SDM using BIOCLIM, Domain, GLM and SVM model profiles

Figure 30. Madii Lii devil’s club Average Score SDM compiled from: BIOCLIM,
Domain, GLM and SVM models.

62

All four models (BIOCLIM, Domain, GLM and SVM; Figure 29), as well as the Average Score
model (Figure 30), suggest that the southern portions of the study area (green coloured area)
provide the most suitable habitat for devil’s club, and much of the northern portions of Madii Lii
predicted to be less suitable or unsuitable habitat. The favourable conditions found in the
southern and western region of the territory corresponds with the lower elevations within the
valleys formed by the various watercourses that flow through the territory and join the Suskwa
River adjacent to the southern boundary.
3.5 Madii lii Predictive Habitat Suitability Generation
3.5.1. Baseline BIOCLIM SDM Model
Though the research had previously created an SDM for existing conditions, to generate a
predictive habitat model encompassing the effects of a warming planet, an alternative baseline
model was constructed using the SSP126:2021–2040 climate model scenario. This BIOCLIM
base model also incorporated the predicted effects of the imminent change in global climate

Figure 31. Madii Lii predicted habitat suitability for devil’s club under current
conditions using BIOCLIM model for the SSP126 2021-2040 climate scenario.

63

temperatures. Comparable to the previous methods the BIOCLIM predictors generated presence
and absence points which then could be used to develop the test and train datapoints (Figure 31).
The BIOCLIM model evaluation (Figure 31) illustrates the training and testing model data
points. Here, the presence (n=18) constitutes the training set which was generated from field data
and are divided into two selections: training and testing. The absence data is pseudo-absence
model evaluation data-points (n=144) data constructed by the software algorithm for measuring
model accuracy. The test-train model received a poor score for the Area Under Curve (AUC)
values of 0.7853 (Table 5).
Table 5. Summary of model evaluation for the BIOCLIM test-train model.
class

Model Evaluation

n presences
n absences
AUC

18
144
0.7853
0.1941
0.0293

correlation coefficient
max TPR+TNR at

3.5.2. MaxEnt Model
The BIOCLIM profile method only considered 'presence' data. The MaxEnt model, on the other
hand, used machine learning to consider data points from presence, absence and background
data. The maximum entropy model predictors were applied to the 17 bioclimatic variables,
which then generated the variable-contribution percentage score for each variable that in turn
were ranked in order of most influence. The results from the contribution plot showed that the
top-ranked environmental variables were Min Temperature of Coldest Month and Mean
Temperature of Coldest Quarter. These two variables provide the most reliable indicator of the
probability of devil’s club presence. These data were then used to generate a series of response
plots for each environmental variable.

64

Figure 32. Multiple response plots for the environmental variables for MaxEnt model generation.

65

These results from the response plots (Figure 32) illustrate the relationship between the
probability of occurrence of devil’s club and each environmental variable. For each plot, the
response is modelled for one environmental variable while the other environmental variables are
held constant at their mean. For example, the largest response is for Bio 13 = Precipitation of
Wettest Month, the probability of devil’s club presence decreases from 70% at 80 mm to 10% at
90 mm precipitation. While for Bio 1 = Annual mean temperature, the response function is a
straight line meaning its effect is negligible and does not affect the model. The results from the
MaxEnt model perform better than the BIOCLIM model, with an AUC of 0.8493 (Table 6). The
AUC of 0.8493 value means there is an 85% chance that the model will reliably distinguish
between positive class and negative classes (Figure 33).

Figure 33. Predictive habitat suitability for devil’s club in Madii Lii for the
SSP245:2021–2040 model using MaxEnt model.

66

Mapped projections of the MaxEnt model output provides a suitable set of model predictions.
Next, the Random Forest model was tested, and output compared with the MaxEnt model to
determine the best model algorithm to use for climate change prediction.
Table 6.Summary of model evaluation scores for MaxEnt model.
class
n presences
n absences
AUC
correlation coefficient
max TPR+TNR at

Model Evaluation
18
144
0.8493
0.4485
0.2341

3.5.3. Random Forest Model
The ROC plot for the Random Forest profile provided an AUC value of 0.9371, (Table 7). In
addition to the MaxEnt and Random Forest, Support Vector Machines models were also applied
to the data. However, the AUC score was only 0.8741, still very good, however for this scenario
Random Forest produced the best performing model and was used to generate the predictive
habitat models for climate change scenarios.
Table 7. Summary of model evaluation scores for Random Forest model.
Class
n presences
n absences
AUC

correlation coefficient
max TPR+TNR at

Model Evaluation
18
144
0.9371
0.6199
0.0612

67

3.6 Predictive Habitat Suitability for Devil’s Club in Madii Lii 2021–2100
The mapped results from the predictive habitat suitability modelling for Madii Lii for SSP245
2021–2100 (moderate/stabilized increase in GHG emissions) are shown below (Figure 34). The
maps depict the expected changes in devil’s club habitat from present-day conditions to the end
of this century in Madii Lii. When comparing 2021–2040 with 2081–2100, it is apparent how
habitat suitability is changing. The plot displaying 2081–2100 has fewer suitable red, yellow, and
green colour. indicating that habitat suitability has declined over the 80-year period. The most
apparent change seen between the 2041–2060 model and the 2061–2080 model. However, there
appears to be an improvement with the amount of suitable habitat in the latter decades of the 21st
century.
When observing the model output for the SSP585:2120–2100 (Figure 35), the red area
(best quality habitat) is considerably reduced when comparing the 2021–2040 map with the
2081–2100 map. From the map it is evident that there is a consistent, gradual decrease in areas of
suitable habitat (red, yellow, and green), which also become more fragmented over time.
Conversely, the areas of poor quality (turquoise) and unsuitable habitat (blue) increase over the
80-year period to the end of the century.

68

Figure 34. Predictive suitability habitat for devil’s club in Madii Lii illustrates the SSP245:2021–
2100 climate.

69

Figure 35. Predictive suitability habitat for devil’s club in Madii Lii illustrates the SSP585:2021–2100 climate.

70

3.6.1. Training and Testing Datapoints with BIOCLIM Model
As well as examining current and future access and availability to devil’s club within Madii Lii,
this research also examined predictive habitat suitability throughout British Columbia. The result
from the BIOCLIM model was used to develop a training and testing presence-absence
algorithm for devil’s club (Figure 36). The model evaluation generated 318 devil’s club present
data-points and 1000 species absent data-points, from the 1600 data-points derived from the
GBIF download for devil’s club locations in BC. The AUC value was 0.8460, providing
satisfactory model to construct the predictive habitat suitability model (Table 8)

Figure 36. British Columbia testing and training points generated to produce BIOCLIM species distribution
model. The blue points represent the species present training data-points and the green the species present
testing data-points, yellow are background training points.

Table 8. Summary of model evaluation scores for BIOCLIM model evaluation.
class
model evaluation
n presences
318
n absences
1000
AUC
0.8460
correlation coefficient
0.4395
.
max TPR+TNR at
0.0228

71

This model is then used to produce the predicted habitat suitability for British Columbia using
the current climate scenario SSP126:2021–2040, (Figure 37).

Figure 37. British Columbia devil’s club suitability habitat using BIOCLIM using current and
predicted near-future conditions for the (SSP126:2021–2040) climate scenario.
3.6.2. Maxent Suitability Model
The BIOCLIM profile method only considers “presence” data. In contrast, the MaxEnt model
profile can generally produce more reliable models using machine learning methods that
consider all data points from presence, absence and background data. The maximum entropy
model was applied to the eighteen environmental variables, which generated predictors that
produced a ranked-percentage variable contribution table. The results from the plot indicate
which environmental variables are the most sensitive to the MaxEnt model. These data showed

72

that the environmental variables, Bio_1 (AMTemp, Annual Mean Temperature) with a 23%
contribution and Bio_2 (MTempD, Mean Temperature of Driest Quarter) with a 22%
contribution, were ranked first and second, respectively, and were the most sensitive to the
model. The results for the variable contribution graph (Appendix 7) and response plot for the
British Columbia Maxent model (Appendix 8). Having determined which contributing
environmental variables influence the predicted likelihood of species presence, the MaxEnt
model was then evaluated for accuracy using the ROC and AUC output. Here the AUC value
produced a value of 0.9176 (Table 9), which is excellent for predicting habitat suitability. When
compared, the BIOCLIM (Figure 38) and the MaxEnt (Figure 39) models showed quite distinct
differences. Overall, the MaxEnt model predicted considerably more suitable habitat across the
region, of which a substantial amount was shown to be good quality devil’s club habitat,
indicated by the red, yellow, and green bands.

Table 9. Summary of MaxEnt model evaluation.
class
model evaluation
n presences
n absences
AUC
correlation coefficient
max TPR+TNR at

318
1000
0.9176
0.6870
0.3892
.

73

Figure 38. British Columbia devil’s club suitability habitat using MaxEnt layer for
SSP126:2021–2040.

3.6.3. Random Forest and SVM Model
The results from the Random Forest and Support Vector Machine model profiles were also tested
and compared using their AUC scores as indictors for best model fit (Appendix 4). The Random
Forest model achieved an AUC score of 0.9720 (Table 10), and SMV received 0.9000.
Therefore, the Random Forest profile was used to complete the predictive habitat suitability
modelling for the SSP245:2021–2100 and SSP585:2021–2100 scenarios. When comparing the
MaxEnt and Random Forest model output maps, the Random Forest model shows smaller

74

scattered areas of high-quality habitat shown in red and suitable areas in yellow, than predicted
by the MaxEnt model.
Table 10. Summary of Random Forest model evaluation.
class
n presences
n absences
AUC

correlation coefficient
max TPR+TNR at

Model Evaluation
318
1000
0.9725
0.8573
0.4173

Figure 39. British Columbia devil’s club suitability habitat using Random Forest layer for
SSP126:2021–2040.

75

3.6.4. Random Forest Predicted Habitat Suitability Model 2021–2100
In the first scenario SSP245:2021-2100, there is only a modest reduction in the amount of the
highest quality habitat shown in red and yellow colour over the four time periods (Figure 40).
The most noticeable change is the decline in the moderate quality habitat (turquoise), which is
replaced with more low-quality habitat (magenta), and there is considerably greater area of nonsuitable habitat (white). Similar to the Madii Lii SSP245 model, the British Columbia PHM
indicates an increase in suitable and highly suitable habitat for devil’s club for the latter part of
the 21st century, possibly due to the plateauing of GHGs in the mid-21st century. SSPP245 states
that there will be an incremental use of green technologies, and therefore a reduction in fossil
fuels emissions. The effect of these improvements will likely take some time to take effect on the
warmer climate, and consequently, devil’s club habitat.
For the second scenario (Figure 41) SSP585:2021–2100, the decline in suitable habitat
for devil's club is severe. In this sequence, there is even less quantity of the highest and very
good habitat (red, yellow, and green) from the outset, and maintains a few very small areas in the
southwest of the province, around the Vancouver area. The loss of mediocre and poorer quality
habitat is quite evident. The scattered area of magenta that runs parallel to the coastline becomes
completely fragmented, leading to infrequent areas of low-quality habitat being replaced by
increasingly large areas of unsuitable habitat as the area of white increases.

76

Figure 40. Predictive habitat suitability models using Random Forest layer for SSP245:2021–2100.

77

Figure 41. Predictive habitat suitability models using Random Forest layer SSP585:2021–2100.

78

Chapter 4
4.0 DISCUSSION
4.1 Looking Back
Madii Lii has been heavily logged for more than sixty years, starting with selective cedar pole
logging in the late 1950s. Since the 1960s, clearcut logging in Madii Lii has removed all
commercial timber from 4169 ha of the landscape, equating to 11.8% of the entire territory—
even more when one accounts for unforested areas in the alpine (Rabnett, 2018). Out of the 352
km2 of territory, most of the 42 km2 of clearcut cutblocks are located within the ICHmc1 and
ICHmc2 BEC zones, which is known to be typical devil’s club habitat. The logging industry has
altered the landscape of Madii Lii, impacting devil’s club habitat and impacting the availability
of large, sought after specimens suitable for medicinal and cultural use.
Research findings suggest that certain types of land use histories significantly affects
devil’s club growth patterns, both in the vigour of individuals and the abundance of devil’s club
stems. Devil’s club stem abundance and stem diameter were both found to be considerably
greater in areas traditionally used and managed by Luutkudziiwus people (i.e., not industrially
logged), compared to both unlogged control sites and previously logged sites. This finding is
significant for Gitxsan harvesters because the most desirable characteristic for devil’s club
harvest is usually its size (stem diameter) and abundance (diversity of choice) as harvest
preference is often based on plant vigour and abundance (Gottesfeld & Anderson, 1988). The
evidence from this study indicates that logging activities are, overall, compromising the
availability and quality of devil's club for Luutkudziiwus people.
Using morphological data as a proxy for devil’s club health and vigour, combined with
stem abundance, the data revealed that devil’s club would take about 26 years to regrow to the
abundance (quantity) and size (diameter) of previously unlogged (control) sites, but at higher

79

elevations, with much slower growing, this may increase to roughly 50 or more years to
regenerate to a size comparable with the preferred individuals harvested for cultural uses, like
those found at the traditional use sites. These data concur with previous research, which suggests
that devil's club could persist after disturbance, such as logging or fire, but regrowth takes a
significant amount of time (Burton & Burton, 2015; Lantz, 2001b; Lantz & Antos, 2002). The
impacts on Luutkudziiwus peoples’ harvesting practices could be significantly impinged with
possibly between 26–50-year recovery period, which is comparable to two generations of people
and therefore affects knowledge transmission. The majority of cultural knowledge regarding
medical plants and healing in Gitxsan contexts is transmitted verbally from one generation to the
next and is practiced experientially on the land. A two-generation gap in knowledge transfer
would likely disrupt the flow of cultural knowledge of traditional medicine (Calixto, 2005).
While ethnobotanical floras and written documents about devil’s club (and other important
plants) are relevant for maintaining knowledge of select plant species, these texts cannot replace
the experiences and relational teachings that come from being on the land and taking part in
harvesting practices directly (e.g., Geniusz, 2022).
Devil’s club growth patterns are not simply a function of the time since logging. As
expected, sites with the most recent logging events had shorter, less productive stems. However,
stem thickness varied with time since industrial logging. For example, morphometric data in
logged sites alone (not in controls or traditionally managed sites) showed that some thicker stems
were observed in sites logged between 1980–1989 and not, as hypothesized, in the older logged
sites from 1970–1979. Devil’s club on sites surveyed from the highest elevations were generally
smaller (containing only a few short, small-diameter stems), making these pockets of devil’s club
especially vulnerable, despite the recovery time. This would suggest that devil’s club grows
particularly slowly at higher elevations. As a result, this phenomenon of smaller devil's club at

80

some of the older sites may be explained by differences in elevation, temperature regime, and
growing season length. However, these data did also show that some of the 1980-89 logging
occurred some of the highest elevations, and that that it is likely there are other factors which
could also influence the stem thickness such as change of logging practices between those eras.
Or it is also possible that devil’s club stem growth is non-linear and may peak at about 40 years,
and possibly decline with stand age
Some practitioners prefer stems harvested from higher elevations, as they are reputedly
more potent or powerful (Larson, 1988; Turner et al., 2011) and “cleaner” compared with plants
harvested from lower elevations (Gottesfeld & Anderson, 1988). Ethnobotanical data suggests
harvesters prefer high elevation sites, possibly because the labour involved offers spiritual
significance or the slower growth could mean denser tissues and higher concentrations of active
ingredients (Armstrong, 2022; Gottesfeld & Anderson, 1988b). This research confirms that even
with a 50+ year recovery period, or even longer at higher altitudes, desired devil’s club harvest
sites are not consistently recovering from logging activities. These findings confirm that forestry
regulators and practitioners need to consider devil’s club recovery on a site-by-site basis and
with regards to cultural preferences. Other research similarly suggests that devil's club is a
negatively impacted understory species that takes longer to recover, regardless of the type of
logging practiced (Beese et al., 2022; Lilles et al., 2018).
During sample collection in Madii Lii, the site elevation and location was noted to be
consistent with the history of logging in the area: intensive high-grading of the most desired
species at lower elevation (easily accessible cedar-hemlock stands) began in the 1950s and
1960s, then the practice of clearcutting rapidly advanced into the territory, removing the most
profitable and accessible wood fibre. By the 1980s, logging occurred at progressively higher

81

elevations (~1000 m), with the continued advancement of cutblocks simultaneously degrading
characteristic devil’s club habitat.
By the later 1990s, loggers in the Kispiox Timber Supply Area (TSA) gradually returned
to the lower elevations and further into the territory to access the more challenging terrain to
harvest in mixed cedar-hemlock forests. This resulted in further destructive impacts on the
remaining devil’s club habitat. As more access roads were constructed, the loggers in the TSA
removed large volumes of timber from more and more parts of the territory, especially at lower
elevations. The rapid construction of access roads and clearcutting meant that most of the easily
accessible, high-economic value forest stands were quickly exhausted, and by the end of the
1990s, 85 % of the wester redcedar stands were removed (Rabnett, 2018).
This research found that the largest and most productive devil’s club patches were
located in the lowest elevation sites. Several factors could result in the larger size and greater
abundance of devil’s club growing at lower elevations. Specimens from lower elevation sites
benefit from longer growing seasons, an earlier snowmelt, and a longer duration of warmer
temperatures compared with those from higher elevation. In addition, it is possible that localized
habitat conditions found at lower elevations such as nitrogen-rich, moist soils with good drainage
all contribute to larger devil’s club patches. The impacts of the early logging should also be
considered. Two different methods of logging, high-grading and clearcutting, both occurred at
these low elevation areas, and may have had different impacts. Those low-elevation logging sites
have had the longest recovery time of 50+ years, possibly affecting the continuity of devil’s club
harvesting and knowledge transmission over recent years.
When preparing devil’s club for medicinal or spiritual purposes, many practitioners
favour the thickest stems for practical reasons: Helen Watkins states “as the labour required to
de-bark and process is the same for a large stem, as it is for a small one” (Betcher, 2016). It

82

makes the task more efficient and produces a greater quantity of medicine for the same amount
of work. The most abundant quantities of thicker stems that can be found in Madii Lii are
generally located at mid- to lower elevations. For example, the site with the thickest diameter
stem (44.9 mm) was collected from an elevation of 471 m within one of the traditional use sites,
among a patch of 29 large diameter stems many of which were 20–25mm diameter. In contrast,
these data showed that at the highest elevation site (at 997 m, logged in 1979) there were only a
handful (13) of stems with a median stem-diameter of 11.4 mm.
Devil’s club showed greater health and vigour at Luutkudziiwus traditional use sites,
such as those close to Gitxsan village sites, the Babine trail network, managed berry sites, and
other cultural sites (e.g., Xsi Nox Nox Loobit). It is likely that these landscapes were managed
and maintained by Luutkudziiwus community members for centuries, resulting in desired
management outcomes for species like devil’s club. Previous research has shown that postlogging slash-burning significantly impaired the regrowth of devil's club (Burton & Burton,
2015; Hamilton, 2006). Factors that impact devil’s club regeneration after slash burning include
altered moisture regime, competition from faster-growing species and, in most cases, highintensity slash burning results in large-scale organic soil depletion. The extremely high
temperatures generated during slash-burning most likely kill all the basal and below-ground buds
and other meristems, severely hindering devil’s club regrowth or survival at these locations.
Other post-logging practices such as chemical treatments (e.g., with the herbicide
glyphosate), and reforestation with quick-growing commercial crops (e.g., spruce and pine), can
also compromise the ability of native species like devil's club to re-establish and persist. Initially,
it was thought that glyphosate’s direct toxic effects were short-lived and limited to the target
species (Grossbard & Atkinson, 1985), however, long-term toxicological studies suggest that

83

plants sprayed with glyphosates experience delayed flowering and disrupted reproduction
patterns (Boutin et al., 2014).
Overall, the logging industry has drastically reduced access to devil’s club in Madii Lii
for Luutkudziiwus community. One of the major factors contributing to peoples’ favourite
patches is accessibility. Geospatial analysis indicates that roughly 150 km2 of Madii Lii territory
lies in the preferred devil’s club habitat of the ICHmc1 and ICHmc2 BEC zones. Of that, ~ 40
km2 have undergone clear-cut logging (~ 27%), and associated impacts from logging road
construction and skidder trails have drastically reduced habitat. Of the remaining 110 km2, 15
km2 of predicted devil’s club habitat is currently inaccessible due to a collapsed bridge that once
spanned Natlan Creek. The collapsed bridge, which would have permitted access to a
considerable area of devil’s club habitat, highlights an important issue that since the forest
licensees ceased operating in Madii Lii, they are no longer responsible for maintaining road
infrastructure. Once roads are no longer needed for harvesting or silviculture, the repair of forest
access roads becomes the duty of the Province, and vehicular access should be provided, where
community stakeholder’s express the desire to maintain access. Lastly, 13 km2 of the remaining
predicted devil’s club habitat, located at the upper reaches of the Denison Creek and the Upper
Suskwa, are also inaccessible due to flooding (likely by a beaver dam).
In addition, the areas of high suitability devil’s club habitat found along the Iltzul Creek
are very steep and are virtually impossible to access. In reality, there are only a few easily
accessible areas of devil's club left in the territory, which are mostly found alongside the Suskwa
Forest Service Road and the Babine Trail. Additional devil’s club patches, which remain
reasonably accessible, are located along the Denison and Iltzul Forest Service Roads. However,
both roads are now heavily congested with brush and are only accessible by all-terrain vehicle or
by foot for those individuals determined to enter these areas. The compounding impacts to

84

habitat, the slow-growing nature of devil's club, cultural preferences, and the territory's
topography mean that access to readily available and desirable stems is considered low.
Results showed that forest canopy and light dynamics affect devil’s club morphology in
significant ways. As predicted, shade had a significant effect on stem length. This makes sense,
as generally plants are light-deprived and stem growth is directed to finding sun flecks and
canopy openings, and therefore, the stems tend to be longer and leggier in areas with less
available sunlight. Partly shaded sites, with mixed canopy gaps, produced devil’s club plants
with large leaves and the thickest diameter stems; the results showed that, where canopy cover
exceeds 75%, many of the smallest diameter stems were located. This also makes sense given
that devil’s club prefers mixed or spotty sun exposure, a characteristic of mature or old growth
forests (MacKinnon & Pojar, 1994). These findings do not fully agree with earlier research
conducted in the Nass Valley, which found that the thickest stems were found in either open
conditions or in the shade created by piles of slash or logs–dead shade (Burton & Burton, 2015).
The difference in findings may be attributed to the differing locations for the research; the Nass
Valley has a more coastal climate with greater cloud, and their sampling occurred in
considerably younger cutblocks where there was more competition with brush or dense trees
(Burton & Burton, 2015).
Previous research performed in second-growth boreal mixed-wood forest showed that
successful understory species such as velvet leaf blueberry (Vaccinium myrtilloides) are able to
persist in heavy shade conditions by employing biomass allocation plasticity. Researchers found
that this species was able to adjust resources to either specific leaf area, leaf area, individual leaf
weight, and the proportion of total biomass in stems, depending on the level of shade (Moola &
Mallik, 1998). Likewise, it appears that devil’s club possesses a similar plasticity trait (found in
most plant species) that induces stem elongation and larger leaf growth under shadier conditions

85

compared with species in open canopy environments (Givnish, 1988; Power et al., 2021;
Valladares & Niinemets, 2008). This research also agrees with previous research which found
devil’s club will flourish in areas exposed to combinations of dense shade and periodic exposure
to bright sunlight, such as found in canopy gaps created by a natural disturbance such as blowdown or lightning strike.
4.2 Looking Forward
The Intergovernmental Panel on Climate Change (IPCC 2018; Masson-Delmotte, et al., 2021)
concluded that anthropogenic climate change has caused an increase in global temperature by
1.18°C above pre-industrial levels and is likely to reach 1.5°C between 2030 and 2050. The
impacts of global temperature increases will vary in different regions and biomes across the
planet. For example, the polar regions are expected to experience increase in average
temperatures of between 5°C to 6°C by 2100 (Price et al., 2013) and the boreal forest biome in
Canada is projected to undergo the widest range of mean annual temperature increases, between
4°C and 11°C by the end of the 21st century (Gauthier et al., 2015; Price et al., 2013). To
understand how these predicted changes may impact devil’s club distribution, range, and habitat
across varying climate scenarios, this research projected species distribution models for varying
climate scenarios with the ultimate goal of predicting how devil’s club will fare under uncertain
climate futures.
Devil’s club SDMs for Madii Lii and British Columbia, combined with climate change
predictions, presented areas of suitable habitat and favourable environmental conditions, which
indicated the probability of devil’s club presence/absence at a particular location. Four high
resolution SDMs in 20-year increments at Madii Lii, and four low resolution SDMs for the
entirety of British Columbia were generated. By extrapolating my documented devil’s club

86

distribution data obtained during this field research SDMs illustrate where devil’s club might
persist (and/or thrive), and where habitat will likely change over the course of the next 80 years.
These data are instrumental for Luutkudziiwus land-use planning, to help identify where devil’s
club habitat should be protected.
Based on these data collected during field sampling in Madii Lii, and the subsequent
PHMs generated from these data-points the climate change scenarios produced some interesting
results. The species distribution models for Madii Lii indicate that steep-sided valleys, toeslopes, and benches, all associated with three large watercourses that flow through Madii Lii
(Natlan, Iltzul and Denison Creeks) and the Suskwa River all continue to provide suitable habitat
for devil's club. The models suggest that two significant indicators for devil’s club
presence/absence are mean annual temperature and mean annual precipitation. The predicted
impacts from both climate change scenarios (SSP245 and SSP585) suggest an overall decline in
highly suitable and suitable habitat over the next 80 years. However there appears to be a distinct
difference in the amount of high-quality habitat at the end of the century between the SSP245
and SSP585 scenarios. The models predict that there will be a considerably larger area of
suitable habitat for devil’s club in the former scenario, most likely due to the lower GHG
emissions as a result of the increased use of green technologies and the decline of fossil fuels
during the mid-21st century. Without those transitions, devil’s club habitat will decrease
significantly in Madii Lii and British Columbia.
The SSP245 model progression shows a decline in habitat through to 2080 but then
begins to rebound, and there is an appreciable increase in the highly-suitable habitat along the
Natlan Creek and Suskwa river areas when comparing 2080 with 2100. The SSP585 model
predicts a steady decline in the areas of highly suitable devil’s club habitat are nearly eliminated,

87

except for a small area by the Suskwa River, and overall, the area of unsuitable habitat increases
over time.
Climate models indicated that annual precipitation and precipitation of the wettest quarter
and driest quarters all had a maximum tolerance, after which precipitation becomes too
excessive, and reduces the probability of devil’s club presence. This will become more
important as precipitation levels become more erratic, and the frequency of extreme-weather
events such as “atmospheric rivers” (ARs) become more common in British Columbia. It is
expected that, between 2070 and 2100, extreme AR events will increase in frequency in British
Columbia, especially during December (Radić et al., 2015). Climate models project an overall
increase in total precipitation over BC, with the increased frequency and intensity of AR-extreme
events, and repeated AR-events with precipitation exceeding 20 mm per day (Radić et al., 2015).
Furthermore, for Madii Lii, the anticipated path for many of these AR-extreme events is to move
northward (from southern BC), bringing more precipitation to the north coast of BC. While
Madii Lii is over 200 km from the coast, there is widespread perception and expectation that the
northwestern region is becoming similar to a coastal in climate, and the additional precipitation
will impact the territory. The northern coastal region of BC is expected to experience annual
precipitation increase between 6% and 26% (dependant on climate change scenario) by the end
of this century (Haughian et al., 2012). Moreover, when comparing the current condition
(SSP126) model with both the intermediate and high-level greenhouse gas emissions scenarios
(SSP 245 and SSP585), my model shows significantly fewer areas of the best and high-value
habitat for devil’s club, with an increasing amount of unsuitable and very low suitability habitats.
However, as with many wild species, there is a progressive northerly shift and expansion of
some previously unsuitable habitats becoming low-quality but suitable (Hirabayashi et al., 2022).
Studies have shown that some species will migrate to higher elevations as the effects of climate

88

change become more apparent (Breshears et al., 2008; Kelly & Goulden, 2008; Lenoir et al.,
2008). Unfortunately, Luutkudziiwus people likely cannot, and will not, have access to these
potential northern clusters at the higher elevations given that they are far from home (reducing
accessibility) or outside of their territory (which would belie customary governance structures
and land-use protocols.). However, it is possible that devil’s club populations from the south and
coastal areas may migrate further inland as the climate changes to a more coastal, as many
researchers are predicting (Haughian et al., 2012). However, with this migration to higher and
mid-elevation habitats may become congested, and competition threatens endemic species with
increased risk of extinction and homogenization (Gómez et al., 2015; Jurasinski & Kreyling,
2007; McKinney & Lockwood, 2001). In all scenarios, the availability and accessibility of
devil’s club in this landscape could be considerably reduced over the next 80 years. Furthermore,
locations where devil’s club clusters persist may be limited to difficult to reach places, such as
remote valley bottoms, and therefore the construction of access roads or trails will be needed to
allow for harvesting. Conversely, any revitalization of traditional trails could benefit
Luutkudziiwus community members as these new trails would possibly serve several purposes
such as resource access and future potential eco-cultural tourism.
Two climate scenarios were analyzed including, SSP 245 and SSP585. Scenario one
(SSP245) first predicts a moderate increase of GHG emissions with resulting warming of 3.84.2oC by the end of the 21st century (Fricko et al., 2017). Scenario two (SSP585) anticipates a
continuous unrestrained worldwide economic growth fuelled by unrestricted burning of hydrocarbons, and thereby producing the highest overall GHG emissions by the year 2100 (Fricko et
al., 2017). As a result, the predicted global temperature increase is expected to be 4.7–5.1oC. by
the end of the 21st century (Fricko et al., 2017). These are graded according to lower impact
(SSP245) and highest impacts (SSP585). The anticipated change in weather patterns could have

89

a profound effect on the access and ubiquity to this vital medicinal and culturally significant
plant for First Nations groups across the province of British Columbia. The baseline model
illustrates the estimated range of devil's club today across BC, and it clearly shows two distinct
regions with good quality habitat, one in the east of the province, on the windward slopes of the
Cariboo, Selkirk, Monashee, and Rocky Mountains, and one in the west, which follows the Coast
Mountains (Figure 42).

Figure 42. British Columbia devil’s club suitability habitat using BIOCLIM using current and
predicted near-future conditions for the (SSP126:2021–2040) climate scenario.

For the predicted models, scenario 1 (SSP245) forecasts a moderate increase in GHGs
leading to ~1.7oC rise in temperatures across BC, while the second scenario (SSP585) predicts a
potential increase of 4.75oC, as GHG emissions continue to rise throughout the 21st century
(Haughian et al. 2012). Climatologists have generated models that calculated the probability of a
90

reoccurrence of the record-breaking temperatures that occurred in British Columbia in June–July
of 2021. The models indicate that by 2055–2064, the chance of an individual day exceeding a
previous record-high is 1-in-1000 years or (0.01%) for the SSP126 (a low emission pathway
reaching net zero by 2075, with global temperatures forecast to increase to by 1.7oC by 2040,
and to 1.3–2.4oC by 2100; Thompson et al., 2022). In contrast, scenario SSP585 projects the
chance of an individual day exceeding a previous recorded high-temperatures (record high)
increasing to 1-in-75 years (1.3%) by 2064, and by 2094, the chance rapidly increases to 1-in-6
years (17%) (Thompson et al., 2022) . This projected increase in temperatures, during the driest
months and throughout the year, could severely impact the viability of a plant that is stressed by
prolonged periods of extreme heat. Combined with the excessive and intense precipitation events
expected, devil’s club is facing a very uncertain future.
Predictions suggest that over the next eighty years, many plant species (not specifically
Oplopanax horridus) will be forced into ever-smaller regions and pushed into areas of higher
elevation and possible higher latitudes (Higgins et al., 2003; Hirabayashi et al., 2022). For most
of BC’s Indigenous peoples, harvesting devil’s club of suitable size and quality will become
increasingly challenging. As mean temperatures increase and weather patterns change, many
understory plant species that flourish in forested microclimates, having cool moist environments
and with lower interannual variability, will experience more pronounced impacts of a warming
climate (De Frenne et al., 2021). Prevéy et al. (2020a) examined three culturally-important foodproducing species from the northwest USA for several Tribes from Washington, Oregon,
northern California, Idaho, and western Montana. These species — beaked hazelnut, Oregon
grape, and salal— are expected to experience a reduction in suitable habitat, affecting cultural
transmission of use and management, and the livelihoods of people that depend on them. The
study also predicts that the projected reduction in climate moisture would cause the southern

91

regions and lower elevation ranges of these species to disappear by the end of the 21st century.
Changes in the thermal regime could also induce a cascade of effects by altering the phenology
of many plant species (Prevéy et al., 2020a).

4.3 Limitations of Research
During this research, I endeavoured to produce a thorough and unbiased assessment of the recent
impacts of logging on devil's club and the predicted impacts of an uncertain climate future.
However, with all research, there are limitations to what was possible due to time and budgetary
restraints. A key factor I would have addressed was that I intended to sample over a wider area.
A significant section of the ICH/clearcut logging areas was inaccessible, and the subsequent
reduction in sampling sites may have skewed my results for both the devil's club abundance and
morphology statistical analyses, as well as the species distribution and climate-change modelling
predictions. Consequently, the results and conclusions derived from this research should be
considered with that caveat in mind. In addition, the Stand-Age/Time Since Disturbance agedata used to date the control and Traditional Use sample sites was retrieved from the Vegetation
Resource Inventory (VRI) (FLNROD, 2018a). Therefore, the reliability of results from the
analysis when compared with the logged sample sites is dependent on the quality of the VRI
data.
Likewise, for the Predicted Habitat Modelling for the Province of British Columbia. The
presence data was retrieved from a single repository (GBIF); therefore, relying on the quality of
that data could be problematic. The GBIF repository depends on the iNaturalist membership to
submit presence data. Consequently, there could be a bias towards recorded presence data points
in areas which have greater populations and accompanying large academic institutions such as

92

universities. Therefore, the areas of highly suitable devil’s club habitat seen in the Vancouver
metropolitan area may be due to a sampling bias.
There are inherent constraints and assumptions involved with predictive habitat
modelling, and model-reliability is derived from the quality of input-data. Models can only
reflect their calibration and training-data; therefore, the entire species-range should be applied to
the model. If a species subset from a limited area is used, then model prediction usually predicts
a contraction of suitable habitat in the future. Another reoccurring issue with PHM is the reliance
on, using the Area Under the receiver operating characteristic Curve (AUC) value as the primary
indicator of model robustness may be problematic, as different spatial scales can result in higher
or lower model fit values (Jiménez-Valverde, 2012; Lobo et al., 2008). However, studies have
found that the Random Forest (RF) model algorithm is often the most effective model and
consistently performs better than the other algorithms at a large spatial scale (Aguirre-Gutiérrez
et al., 2013).
PHMs cannot discriminate the effects of inter-species influences and do not incorporate
many other environmental factors such as soil nutrient composition, mycorrhizal relationships,
and soil moisture, which could induce either a positive or negative effect on a species
distributions and adaptations to climate change. Moreover, SDMs and PHMs also typically fail
to take account of population differences in their genetic adaptation to local climates, soils, and
biotic stresses. In addition, PHMs do not consider the impact of change of land-use or
disturbances such as fire or beetle infestation. Refugia from such disturbances, and locations
more likely to provide habitat continuity over time can also offer protection for suitable habitats
that otherwise might be lost (Rose & Burton, 2009). Despite these limitations, SDMs and PHMs
modelling provide valuable insight into our likely near-future reality and provides a basis for
scenario planning and can guide the direction of conservation efforts.

93

Chapter 5
5.0 CONCLUSION
The primary goal of this master’s research was to assess the effects of a history of clearcutting on
the wa’uumst, or devil’s club (Oplopanax horridus) resource in the wilp Luutkudziiwus
traditional territory known as Madii Lii, and to understand how this impacts the current and
future availability of this medicinal plant. I examined whether devil's club can persist and regrow
after the logging, and what the effects are on plant abundance and morphology in former
clearcuts compared to undisturbed and/or Indigenous managed areas. Predictive habitat models
using two climate change scenarios (SSP245 and SSP585) were used to assess local and regional
climate change on devil’s club distribution in British Columbia.
The results from this research suggest that devil's club has been significantly impacted by
clearcut logging practices in Madii Lii over the last five decades, with reduced abundance, size
and availability. However, as with previous studies (Burton & Burton, 2015; Lantz and Antos,
2002), this research also indicates that devil's club can persist in cutblocks after logging.
Unfortunately, it appears that to reach the preferred sizes (e.g., >25 mm stem diameter, Burton &
Burton, 2015) and quantities for medicinal and cultural purposes may take 26 to possibly 50 or
more years to recover, depending on elevation. The research also showed that the optimum shade
where devil’s club can thrive is up to 75%, suggesting post-logging recovery could be impeded
where replanted with closely packed commercial crops such as hybrid-spruce replacing the
original cedar-hemlock dominated forests.

94

5.1 Recommendations
Given the cultural significance and the potential impact on Luutkudziiwus and the loss of
valuable medicinal knowledge, any future extractive or destructive resource projects should
include protective measures for areas where devil's club is regularly harvested and patches of
devil’s club which exceed a pre-determined size. To protect these valuable patches, I suggest
determining a buffer zone of a proportionate size using green-tree retention areas (Bennett &
Mulongoy, 2006; Hansson, 2001; Rosenvald & Lõhmus, 2008). Madii Lii is currently under
threat from linear developments such as liquified natural gas (LNG) pipelines and power
transmission lines, including the Prince Rupert Gas Transmission project (PRGT) and, more
recently, the possible revival of PRGT to service the proposed Ksi Lisims LNG project. If these
or similar developments eventually proceed, then it is essential that exclusion zones surrounding
productive devil’s club habitat should be implemented prior to commencement of any work.
Another concern for continued cultural use is that only a handful of the areas sampled
had a plentiful number of stems of the favoured diameter and are possibly vulnerable to overharvest. Because of the slow-growing nature of devil’s club, and the prolonged time required to
reach the desired size, the low availability of high-quality devil’s club stems may inhibit the flow
of cultural knowledge of devil’s club. Therefore, by utilizing these data from survey plots and
SDMs for devil’s club in Madii Lii, it would be possible for Luutkudziiwus Elders or knowledge
keepers to regulate the quantity and location of any significantly large devil’s club harvest in
Madii Lii. Also, in order to prevent a catastrophic loss of adequately sized devil’s club, it may
be wise to prepare for future scarcity by the initiation of a devil’s club cultivation plan. One
possibility is to develop a multi-pronged devil’s club sustainability plan, based on protection,
regulated harvesting, and possibly cultivation to be implemented by Luutkudziiwus community

95

leaders. Cultivation of devil’s club could occur within a forest garden, or in areas specifically
designated within Madii Lii. such as proximate to the ancient village site. This may not be the
ideal solution, due to the spiritual nature of harvesting but may ultimately be required to ensure
the continued availability of this important medicinal plant.
The historic victory in the Federal Supreme Court Decision of Delgamuukw and
Gisday'wa v. the Crown (S.C.C.44) affirmed Aboriginal title in British Columbia, that it is a right
to the land itself, and not only the right to hunt, fish or harvest. The ruling further means that
when dealing with Crown land, the government must consult with and may have to compensate
title holders whose rights are affected. The duty to consult with First Nations groups affected by
natural resource projects, including forestry and pipelines, is often seen as a hurdle by industry
rather than an opportunity to understand how these projects and industrial activities will likely
impact traditional livelihoods including the gathering of plants for food and medicine, in this
case, devil’s club. As a plaintiff in the case, Luutkudziiwus has not been effectively consulted in
resource management contexts, and this research confirms how impacts from industry directly
impact the cultural and ecological health and well-being of house members.
As well as the data regarding the location and size of devil’s club patches, this research
complements the catalogue of evidence that details how successive provincial governments and
their agents have largely ignored the incidental impacts of intensive forest practices such as
widespread clearcutting. The detrimental impact of even-aged forest management on devil's club
and these effects on Luutkudziiwus members, have been effectively disregarded by forest
licensees and the provincial government for the last fifty years. Fortunately, as this research
shows, devil's club can recover from the impacts of logging, though it appears that several
decades, or two generations of Luutkudziiwus will pass during this time.

96

In conclusion, possibly the most significant contribution of this study is the
documentation of the detrimental effects of environmental change on the culturally salient
medicinal plant known as wa’uumst, or devil’s club in Madii Lii. As the highly suitable devil’s
club habitat degrades and becomes less available due to the warming climate, devil's club of
appropriate size will likely become more challenging to find and harvest, particularly with
respect to the SSP585 climate model. Consequently, the flow of cultural knowledge may
diminish with the availability of the plant, both in Madii Lii and throughout British Columbia.
This environmental degradation has seen immediate impacts in the case of logging, and this
degradation may continue over time as a consequence of climate change. The reduction in
availability and the impaired access to devil's club has and will continue to harm the transfer of
Luutkudziiwus Indigenous and local knowledge.

5.2 Future Research
These results illustrate the time required for devil’s club to re-establish post-logging, and,
therefore, how this may impact the availability of the plant for Luutkudziiwus cultural use.
However, as time, resources and access were limiting factors, for future research I suggest
undertaking a greater research area throughout the numerous Gitxsan wilp territories and
biogeoclimatic site series locations. As for the question regarding shade, there appears to be
varied result regarding the impact of canopy openness on the morphology and vigour of devil's
club. Therefore, I suggest a long-term controlled experiment (e.g., 5–10 years), where shade
screens control different amounts of shade and the plants' response to the treatment is measured.
A similar controlled experimental technique could be used to measure the response for different
soil moisture contents and different soil nutrient levels (e.g., nitrogen).

97

During this research I also examined the relationships between devil's club and the
human activities, especially logging, on Luutkudziiwus traditional territory of Madii Lii. The
interpretation of the results suggests that devil’s club can persist and re-establish post-logging.
However, as previously stated, the time required for this process to occur is still variable, which
may be a result of influences in elevation and other site attributes. Therefore, a wide-ranging
study of different site types from different regions will help determine if there is a definitive
time-scale for devil’s club recovery after clearcut logging to occur. Possible important factors
that have been examined are soil moisture and soil nutrients. I intended to measure the influence
of soil moisture on devil’s club fecundity, but the soil-moisture sensing equipment failed and this
part of the research was made impossible. Therefore, any future research should evaluate the role
of soil moisture and soil nutrient levels in the productivity of devil’s club within a logged
environment. In addition, I suggest that future research should be expanded to incorporate the
use and monitoring of the effects of partial cutting and gap-based silviculture, as employed in the
Date Creek Experimental Forest near Kispiox, and as described in Lilles et al. (2018).

98

Literature Cited

Aguirre-Gutiérrez, J., Carvalheiro, L. G., Polce, C., Loon, E. E. van, Raes, N., Reemer, M., &
Biesmeijer, J. C. (2013). Fit-for-purpose: Species distribution model performance
depends on evaluation criteria – dutch hoverflies as a case study. PLOS ONE, 8(5),
e63708. https://doi.org/10.1371/journal.pone.0063708
Alaback, P. B. (1984). Plant succession following logging in the Sitka spruce-western hemlock
forests of southeast Alaska: Implications for management (Vol. 173). US Department of
Agriculture, Forest Service, Pacific Northwest Forest and Range Experiment Station.
(Vol. 173).
Applequist, W. L., Brinckmann, J. A., Cunningham, A. B., Hart, R. E., Heinrich, M., Katerere,
D. R., & Van Andel, T. (2020). Scientistsʼ warning on climate change and medicinal
plants. Planta Medica, 86(01), 10–18.
Armstrong, C. G. D. (2017). Hist orical ecology of cultural landscapes in the Pacific Northwest
[PhD Thesis]. Environment: Department of Archaeology. Simon Fraser University.
Burnaby, British Columbia.
Armstrong CG. 2022. Silm Da’axk (to Revive and Heal Again): Historical Ecology and
Ethnobotany in Gitselasu Lahkhyuup. Mitchell Press, Vancouver, BC.
Banner, A., (1993). A field guide to site identification and interpretation for the Prince Rupert
Forest Region (Vol. 26). Ministry of Forests, Research Program..
Barker, J. S., Simard, S. W., & Jones, M. D. (2014). Clearcutting and high severity wildfire have
comparable effects on growth of direct-seeded interior Douglas-fir. Forest Ecology and
Management, 331, 188–195. https://doi.org/10.1016/j.foreco.2014.08.004

99

BC Ministry of Forests and Range. (2009). Administrative Guide for Forest Stewardship Plans:
Volume 1: Preparation and Approval of an FSP Version 2.1. BC Ministry of Forests and
Range.
Beese, W. J., Sandford, J. S., Harrison, M. L., & Filipescu, C. N. (2022). Understory vegetation
response to alternative silvicultural systems in coastal British Columbia montane forests.
Forest Ecology and Management, 504, 119817.
https://doi.org/10.1016/j.foreco.2021.119817
Benner, J., Knudby, A., Nielsen, J., Krawchuk, M., & Lertzman, K. (2019). Combining data
from field surveys and archaeological records to predict the distribution of culturally
important trees. Diversity and Distributions, 25(9), 1375–1387.
https://doi.org/10.1111/ddi.12947
Bennett, G., & Mulongoy, K. J. (2006). Review of experience with ecological networks,
corridors and buffer zones. Secretariat of the Convention on Biological Diversity,
Montreal, Technical Series, 23, 100.
Betcher. (2016). Devils Club; Tlingit Traditions of Helen Watkins. Ties To Alaska's Wild Plants
series [Film]. University of Alaska. Museum of the North. Fairbanks. Alaska. USA.
[Video]https://www.youtube.com/watch?v=DyUssKc2TLQhttps://www.youtube.com/wa
tch?v=DyUssKc2TLQ

Blanquet, B. (1964). Grudzuge der Vegetationskunde. Springer Verlag, Wien.
Blaschke. (1842). Topographia Medica Portus Novi-Archangelscensis. Wiehoberi, St.
Petersburg.

100

Bogawski, P., Damen, T., Nowak, M. M., Pędziwiatr, K., Wilkin, P., Mwachala, G.,
Pierzchalska, J., & Wiland-Szymańska, J. (2019). Current and future potential
distributions of three Dracaena Vand. Ex L. species under two contrasting climate change
scenarios in Africa. Ecology and Evolution, 9(12), 6833–6848.
Boutin, C., Strandberg, B., Carpenter, D., Mathiassen, S. K., & Thomas, P. J. (2014). Herbicide
impact on non-target plant reproduction: What are the toxicological and ecological
implications? Environmental Pollution, 185, 295–306.
Bowman, D. M., & Haberle, S. G. (2010). Paradise burnt: How colonizing humans transform
landscapes with fire. Proceedings of the National Academy of Sciences, 107(50), 21234–
21235.
Boyle, C. A., Lavkulich, L., Schreier, H., & Kiss, E. (1997). Changes in land cover and
subsequent effects on Lower Fraser Basin ecosystems from 1827 to 1990. Environmental
Management, 21(2), 185–196.
Breiman, L. (2001). Random forests. Machine Learning, 45(1), 5–32.
Breshears, D. D., Huxman, T. E., Adams, H. D., Zou, C. B., & Davison, J. E. (2008). Vegetation
synchronously leans upslope as climate warms. Proceedings of the National Academy of
Sciences, 105(33), 11591–11592.
British Columbia Ministry of Forests, Lands, Natural Resource Operations and Rural
Development. 2018a. Forest analysis and inventory: VRI — forest vegetation composite
polygons and layer 1. Available from https://catalogue. data.gov.bc.ca/dataset/vri-2019forest-vegetation-composite-polygons [accessed September 18th 2020].
British Columbia Ministry of Forests, Lands, Natural Resource Operations and Rural
Development. 2018b. Forest analysis and inventory: harvested areas of BC (consolidated

101

cutblocks). Available from https://catalogue.data.gov.bc.ca/ dataset/harvested-areas-ofbc-consolidated-cutblocks- [accessed 18 September 2020].

British Columbia Ministry of Forests, Lands, Natural Resource Operations and Rural
Development. 2017a. Forest tenures. FADM: Timber Supply Area (TSA). Available
from https://catalogue.data.gov.bc.ca/dataset/fadm-timber-supply- area-tsa [accessed
September 20th 2020].

British Columbia Ministry of Forests, Lands, Natural Resource Operations and Rural
Development. 2018c. Forest tenures: forest tenure road section lines. Available from
https://catalogue.data.gov.bc.ca/dataset/forest-tenure-road- section-lines [accessed 18
September 2020].
British Columbia Ministry of Forests, Lands, Natural Resource Operations and Rural
Development. 2018d. GeoBC: Freshwater Atlas Stream Network. Available from
https://catalogue.data.gov.bc.ca/dataset/freshwater-atlas-stream- network [accessed
September 20th 2020].
British Columbia Ministry of Forests and Range. 2021. Research Branch. Biogeoclimatic
Ecosystem Classification Program. Province of British Columbia, Victoria, B.C.
Available from https://www.for.gov.bc.ca/hre/becweb/index.html. [Accessed 15
September 2021]
British Columbia Species & Ecosystems Explorer (2023). BC Conservation Data Centre.
Province of British Columbia, Victoria, B.C. Available from:https://a100.gov.bc.ca/pub/eswp/search.do [accessed April 28th 2023].

102

British Columbia. Ministry of Water, Land and Air Protection. 2004. Western redcedar–
Douglas-fir/devil’s-club (Thuja plicata-Pseudotsuga menziesii-Oplopanax horridus) in
Accounts and Measures for Managing Identified Wildlife – Accounts V. 2004. B.C.
Ministry of Water, Land and Air Protection, Victoria, B.C. Available from:
http://www.env.gov.bc.ca/wld/frpa/iwms/accounts.html accessed April 28th 2023]
Buenz, E. J. (2005). Country development does not presuppose the loss of forest resources for
traditional medicine use. Journal of Ethnopharmacology, 100(1), 118–123.
https://doi.org/10.1016/j.jep.2005.05.005
Bunce, R. G. H., Bogers, M. M. B., Evans, D., Halada, L., Jongman, R. H. G., Mucher, C. A.,
Bauch, B., De Blust, G., Parr, T. W., & Olsvig-Whittaker, L. (2013). The significance of
habitats as indicators of biodiversity and their links to species. Ecological Indicators, 33,
19–25.
Burton, C. M., & Burton, P. J. (2015). Recovery of Oplopanax horridus (Sm.) Miq., an Important
Ethnobotanical Resource, after Clearcut Logging in Northwestern British Columbia.
Ethnobotany Research and Applications, 14, 001. https://doi.org/10.17348/era.14.0.001015
Burton, P., Burton, C., & McCulloch, L. (2000). Exploring options for the management of wild
berries in the Kispiox Forest District: Phase one of a pilot project focusing on the Suskwa
River area. Prepared for Min. Forests, Smithers, BC. Retrieved on the Web at
Http://Sernbc. ca/Pdf/Berry_Management/Options% 20for% 20Management% 20of%
20Wild% 20Ber Ries% 20in% 20KFD% 20Final% 20Report. Pdf.
Calixto, J. B. (2005). Twenty-five years of research on medicinal plants in Latin America: A
personal view. Journal of Ethnopharmacology, 100(1–2), 131–134.

103

Calway, T., Du, G.-J., Wang, C.-Z., Huang, W.-H., Zhao, J., Li, S.-P., & Yuan, C.-S. (2012).
Chemical and pharmacological studies of Oplopanax horridus, a North American
botanical. Journal of Natural Medicines, 66(2), 249–256.
Canada Natural Resources. (2019). High Resolution Digital Elevation Model (HRDEM)–
CanElevation Series–Product specifications Government of Canada Natural Resources
Canada Ottawa. https://open.canada.ca/en/open-maps
Carrero, G. C., Fearnside, P. M., do Valle, D. R., & de Souza Alves, C. (2020). Deforestation
trajectories on a development frontier in the Brazilian Amazon: 35 years of settlement
colonization, policy and economic shifts, and land accumulation. Environmental
Management, 66(6), 966–984.
Chandler, J. R., Haeussler, S., Hamilton, E. H., Feller, M., Bradfield, G., & Simard, S. W.
(2017). Twenty years of ecosystem response after clearcutting and slashburning in
conifer forests of central British Columbia, Canada. PLOS ONE, 12(2), e0172667.
https://doi.org/10.1371/journal.pone.0172667
Cheung, S. S., Hasman, D., Khelifi, D., Tai, J., Smith, R. W., & Warnock, G. L. (2019). Devil’s
Club falcarinol-type polyacetylenes inhibit pancreatic cancer cell proliferation. Nutrition
and Cancer, 71(2), 301–311.
Cheung, S. S., Tai, J., Hasman, D., Ou, D., & Warnock, G. L. (2015). Inhibition of human
pancreatic cancer cell proliferation by Devil’s Club Oplopanax horridus and its
Polyacetylene bioactive compound. Nutrition and Cancer, 67(6), 954–964.
Coates, K. D., & Burton, P. J. (1997). A gap-based approach for development of silvicultural
systems to address ecosystem management objectives. Forest Ecology and Management,
99(3), 337–354

104

Compton. D. B. (1993). Upper North Wakashan and Southern Tsimshian ethnobotany: The
knowledge and usage of plants and fungi among the Oweekeno, Hanaksiala (Kitlope and
Kemano), Haisla (Kitamaat) and Kitasoo Peoples of the Central and North Coasts of
British Columbia Doctoral dissertation, University of British Columbia Vancouver BC.
Condro, A. A., Prasetyo, L. B., Rushayati, S. B., Santikayasa, I. P., & Iskandar, E. (2021).
Predicting hotspots and prioritizing protected areas for endangered primate species in
Indonesia under changing climate. Biology, 10(2), 154.
Cruse‐Sanders, J. M., & Hamrick, J. L. (2004). Genetic diversity in harvested and protected
populations of wild American ginseng, Panax quinquefolius L.(Araliaceae). American
Journal of Botany, 91(4), 540-548
Curtis E. The Haida. In: The North American Indian. New York: Johnson Reprint, 1916: 115210. (S. Hodge., Ed. The Indians of the United States, the Dominion of Canada, and
Alaska; (Vol. 2).
Dawe, C. A., Filicetti, A. T., & Nielsen, S. E. (2017). Effects of linear disturbances and fire
severity on velvet leaf blueberry abundance, vigor, and berry production in recently
burned jack pine forests. Forests, 8(10), 398.
Dawson, N., Coolsaet, B., Sterling, E., Loveridge, R., Gross-Camp, Nicole D., Wongbusarakum,
S., Sangha, K., Scherl, L., Phan, H. P., Zafra-Calvo, N., Lavey, W., Byakagaba, P.,
Idrobo, C. J., Chenet, A., Bennett, N., Mansourian, S., & Rosado-May, F. (2021). The
role of Indigenous peoples and local communities in effective and equitable conservation.
Ecology and Society, 26(3). https://doi.org/10.5751/ES-12625-260319
De Frenne, P., Lenoir, J., Luoto, M., Scheffers, B. R., Zellweger, F., Aalto, J., Ashcroft, M. B.,
Christiansen, D. M., Decocq, G., De Pauw, K., Govaert, S., Greiser, C., Gril, E., Hampe,

105

A., Jucker, T., Klinges, D. H., Koelemeijer, I. A., Lembrechts, J. J., Marrec, R., …
Hylander, K. (2021). Forest microclimates and climate change: Importance, drivers and
future research agenda. Global Change Biology, 27(11), 2279–2297.
https://doi.org/10.1111/gcb.15569
de Groot, A. (2005). Review of the Hydrology, Geomorphology, Ecology and Management of the
Skeena River Floodplain. (Drosera Ecological Consulting) [Bulkley Valley Centre for
Natural Resources Research & Management] Smithers BC.
https://bvcentre.ca/files/research_reports/04-03SkeenaIslandsReview.pdf
Deagle, G. (1988). Traditional West Coast Native Medicine. Canadian Family Physician, 34,
1577–1580.
Delgamuukw v. British Columbia File No 23799. 1997. Reasons for Judgment of the Supreme
Court of Canada [also known as Delgamuukw III]
Dormann, C. F. (2007). Promising the future? Global change projections of species distributions.
Basic and Applied Ecology, 8(5), 387–397.
Drummond, C. P. (2018). Phylogeography of the Plants Disjunct between Western North
America and the Great Lakes Region: Patterns, Case Studies, and Implications for
Future Comparative Research. The University of Wisconsin-Madison.
Elahmer, N. (2018). Bioassay-guided antidiabetic potentials of Devil’s club (Oplopanax
horridus) preparations from the traditional pharmacopeia of the Squamish and other first
nations of British Columbia [Masters (MSc), University of Montreal].
https://papyrus.bib.umontreal.ca/xmlui/bitstream/handle/1866/21369/Nyruz_Elahmer_20
18_Memoire.pdf?sequence=2&isAllowed=y

106

ESRI. (2020). ArcGIS Desktop: Redlands, CA: Environmental Systems Research Institute.
(Release 10.8.1). ESRI.
Fa, J. E., Watson, J. E., Leiper, I., Potapov, P., Evans, T. D., Burgess, N. D., Molnár, Z.,
Fernández-Llamazares, Á., Duncan, T., & Wang, S. (2020). Importance of Indigenous
Peoples’ lands for the conservation of Intact Forest Landscapes. Frontiers in Ecology and
the Environment, 18(3), 135–140.
Farnsworth, N. R., & Soejarto, D. D. (1991). Global importance of medicinal plants. The
Conservation of Medicinal Plants, 26, 25–51.
Fick, S. E., & Hijmans, R. J. (2017). WorldClim 2: New 1‐km spatial resolution climate surfaces
for global land areas. International Journal of Climatology 37:4302-4315.
www.worldclim.org .
Franklin, J. (2010). Mapping species distributions: Spatial inference and prediction. Cambridge
University Press.
Fricko, O., Havlik, P., Rogelj, J., Klimont, Z., Gusti, M., Johnson, N., Kolp, P., Strubegger, M.,
Valin, H., & Amann, M. (2017). The marker quantification of the Shared Socioeconomic
Pathway 2: A middle-of-the-road scenario for the 21st century. Global Environmental
Change, 42, 251–267.
Garibaldi, A., & Turner, N. (2004). Cultural keystone species: Implications for ecological
conservation and restoration. Ecology and Society, 9(3).
Garnett, S. T., Burgess, N. D., Fa, J. E., Fernández-Llamazares, Á., Molnár, Z., Robinson, C. J.,
Watson, J. E., Zander, K. K., Austin, B., & Brondizio, E. S. (2018). A spatial overview of
the global importance of Indigenous lands for conservation. Nature Sustainability, 1(7),
369–374.

107

Gauthier, S., Bernier, P., Kuuluvainen, T., Shvidenko, A. Z., & Schepaschenko, D. G. (2015).
Boreal forest health and global change. Science, 349(6250), 819–822.
.Global Biodiversity Information Facility.(2022.) [GBIF.org]. www.gbif.org
Gedney, N., & Valdes, P. J. (2000). The effect of Amazonian deforestation on the northern
hemisphere circulation and climate. Geophysical Research Letters, 27(19), 3053–3056.
Geertsema, M., Cruden, D. M., & Schwab, J. W. (2006). A large rapid landslide in sensitive
glaciomarine sediments at Mink Creek, northwestern British Columbia, Canada.
Engineering Geology, 83(1–3), 36–63.
Geniusz, W. M. (2022). Our knowledge is not primitive: Decolonizing botanical Anishinaabe
teachings. Syracuse University Press.
Getman-Pickering, Z. L., Campbell, A., Aflitto, N., Grele, A., Davis, J. K., & Ugine, T. A.
(2020). LeafByte: A mobile application that measures leaf area and herbivory quickly
and accurately. Methods in Ecology and Evolution, 11(2), 215–221.
Givnish, T. J. (1988). Adaptation to sun and shade: A whole-plant perspective. Functional Plant
Biology, 15(2), 63–92. https://doi.org/10.1071/pp9880063
Gómez, J. M., González-Megías, A., Lorite, J., Abdelaziz, M., & Perfectti, F. (2015). The silent
extinction: Climate change and the potential hybridization-mediated extinction of
endemic high-mountain plants. Biodiversity and Conservation, 24(8), 1843–1857.
https://doi.org/10.1007/s10531-015-0909-5
Gordon, D., Lorenz, A., & Friesen, M. (1996). Lakelse WRP Project, Level 1 fisheries
assessment. Terrace, BC.

108

Gottesfeld, A. S., Rabnett, K. A., & Hall, P. E. (2002). Conserving Skeena fish populations and
their habitat. Skeena Fisheries Commission, Box, Box 229, Hazelton BC.
https://psf.ca/wp-content/uploads/2021/10/Download-PDF561-1.pdf
Gottesfeld, L. M. J., & Anderson, B. (1988a). Gitksan traditional medicine: Herbs and healing.
Journal of Ethnobiology, 8(1), 13–33.
Gregory, S. D., Brook, B. W., Goossens, B., Ancrenaz, M., Alfred, R., Ambu, L. N., & Fordham,
D. A. (2012). Long-term field data and climate-habitat models Show that orangutan
persistence depends on effective forest management and greenhouse gas mitigation.
PLOS ONE, 7(9), e43846. https://doi.org/10.1371/journal.pone.0043846
Grossbard, E., & Atkinson, D. (1985). The herbicide glyphosate (Vol. 2). Butterworths London.
Guisan, A., & Zimmermann, N. E. (2000). Predictive habitat distribution models in ecology.
Ecological Modelling, 135(2–3), 147–186.
Haeussler, S., Canham, C. D., & Coates, K. D. (2013). Complexity in temperate forest dynamics.
(pp. 74–92) in C. Messier, K. J. Puettmann & K. D. Coates (Eds.). Managing forests as
complex adaptive systems: building resilience to the challenge of global change.
Routledge
Haeussler, S. & D. Kneeshaw, 2003. Comparing forest management to natural processes. pages
307-368 in P. J. Burton, C. Messier, D. W. Smith & W. L. Adamowicz (eds.). Towards
Sustainable Management of the Boreal Forest. NRC Research Press, Ottawa, Ontario.
Hamilton, E. H. (2006). Fire effects and post-burn vegetation development in the Sub-Boreal
Spruce zone: Mackenzie (Windy Point) site. BC Min. For. Range, Res. Br., Victoria.
Hansson, L. (2001). Key habitats in Swedish managed forests. Scandinavian Journal of Forest
Research, 16(sup003), 52–61. https://doi.org/10.1080/028275801300090609

109

Harper, K. A., Macdonald, S. E., Burton, P. J., Chen, J., Brosofske, K. D., Saunders, S. C.,
Euskirchen, E. S., Roberts, D. A. R., Jaiteh, M. S., & Esseen, P.-A. (2005). Edge
influence on forest structure and composition in fragmented landscapes. Conservation
Biology, 19(3), 768–782.
Haughian, S. R., Burton, P. J., Taylor, S. W., & Curry, C. (2012). Expected effects of climate
change on forest disturbance regimes in British Columbia. Journal of Ecosystems and
Management, 13(1).
Hemstrom, M. A., & Logan, S. E. (1986). Plant association and management guide: Siuslaw
National Forest (Gen. Tech. Rep. PNW-GTR-R6-Eco 220.). US Department of
Agriculture, Forest Service, Pacific Northwest Region.
Heritage Conservation Act. RSBC 1996, c187.
Higgins, S. I., Lavorel, S., & Revilla, E. (2003). Estimating plant migration rates under habitat
loss and fragmentation. Oikos, 101(2), 354–366.
Hijmans, R. J., & Elith, J. (2021). Species distribution models. https://rspatial. org/sdm
Hirabayashi, K., Murch, S. J., & Erland, L. A. (2022). Predicted impacts of climate change on
wild and commercial berry habitats will have food security, conservation and agricultural
implications. Science of The Total Environment, 157341.
Hoffman, K. M., Davis, E. L., Wickham, S. B., Schang, K., Johnson, A., Larking, T., Lauriault,
P. N., Quynh Le, N., Swerdfager, E., & Trant, A. J. (2021). Conservation of Earth’s
biodiversity is embedded in Indigenous fire stewardship. Proceedings of the National
Academy of Sciences, 118(32), e2105073118.
Iverson, L. R., & Prasad, A. M. (2001). Potential changes in tree species richness and forest
community types following climate change. Ecosystems, 4(3), 186–199.

110

Jiménez-Valverde, A. (2012). Insights into the area under the receiver operating characteristic
curve (AUC) as a discrimination measure in species distribution modelling. Global
Ecology and Biogeography, 21(4), 498–507. https://doi.org/10.1111/j.14668238.2011.00683.x
Johnson, A. (2020). Medicinal, cultural, and spiritual relationships between British Columbia
First Nations and Oplopanax horridus. SLC Undergraduate Writing Contest, 3.
Johnson, L. M. (2006a). Gitksan medicinal plants-cultural choice and efficacy. Journal of
Ethnobiology and Ethnomedicine, 2(1), 1–23.
Johnson, L. M. (2006b). Gitksan medicinal plants-cultural choice and efficacy. Journal of
Ethnobiology and Ethnomedicine, 2(1), 1–23.
Johnson, L. M. (2013). Plants, places, and the storied landscape: Looking at First Nations
perspectives on plants and land. BC Studies: The British Columbian Quarterly, 179, 85–
105.
Jurasinski, G., & Kreyling, J. (2007). Upward shift of alpine plants increases floristic similarity
of mountain summits. Journal of Vegetation Science, 18(5), 711–718.
Kachel, S., Herbert, J., & Banner, C. (2006). Ethnobotany and native plant production.
Ethnobotanically significant plants of the Pacific. University of Washington, Seattle.
https://depts.washington.edu/propplnt/Chapters/Ethnobotany%20&%20Native%20Plant
%20Production.pdf
Kaky, E., & Gilbert, F. (2016). Using species distribution models to assess the importance of
Egypt’s protected areas for the conservation of medicinal plants. Journal of Arid
Environments, 135, 140–146. https://doi.org/10.1016/j.jaridenv.2016.09.001

111

Kaye, D. (2012). Prince George Defined Forest Area Sustainable Forest Management Plan.
Canadian Forest Products Ltd. Vancouver, British Columbia.
Keenan, R. J., & Kimmins, J. P. (1993). The ecological effects of clear-cutting. Environmental
Reviews, 1(2), 121–144.
Kelly, A. E., & Goulden, M. L. (2008). Rapid shifts in plant distribution with recent climate
change. Proceedings of the National Academy of Sciences, 105(33), 11823–11826.
Klinkenberg, B. (2021). E-Flora BC: Electronic Atlas of the Flora of British Columbia
[eflora.bc.ca]. Lab for Advanced Spatial Analysis, Department of Geography, University
of British Columbia, Vancouver. [Accessed Dec 12th 2021].
Komsta, L. (2013). mblm: Median-based linear models. R Package Version 0.12, 1.
Kramkowski, V. (2012). END-games: The material and political interests of emulating natural
disturbances in the Canadian boreal forest. Environments, 38(1), 81–98.
Lantz, T. C. (2001a). Examining the potential role of co-operatives in the ethical
commercialization of medicinal plants: plant conservation, intellectual property rights,
ethics, and devil's club (Oplopanax horridus). BC Institute for Co-operative Studies,
University of Victoria. British Columbia

Lantz, T. C. (2001). The Population Ecology and Ethnobotany of Devil's Club (Oplopanax
Horridus (Sm.) Torr. & A. Gray Ex. Miq.; Araliaceae) (Doctoral dissertation, University
of Victoria)
Lantz, T. C., & Antos, J. A. (2002). Clonal expansion in the deciduous understory shrub, devil’s
club (Oplopanax horridus; Araliaceae). Canadian Journal of Botany, 80(10), 1052–1062.

112

Lantz, T. C., Swerhun, K., & Turner, N. J. (2004). Devil’s club (Oplopanax horridus): An
ethnobotanical review. HerbalGram.
Larson, R. A. (1988). The antioxidants of higher plants. Phytochemistry, 27(4), 969–978.
Lawler, J. J., Shafer, S. L., White, D., Kareiva, P., Maurer, E. P., Blaustein, A. R., & Bartlein, P.
J. (2009). Projected climate-induced faunal change in the Western Hemisphere. Ecology,
90(3), 588–597.
Lenoir, J., Gégout, J.-C., Marquet, P. A., de Ruffray, P., & Brisse, H. (2008). A significant
upward shift in plant species optimum elevation during the 20th century. Science,
320(5884), 1768–1771.
Lepofsky, D., Armstrong, C. G., Mathews, D., & Greening, S. (2020). Understanding the past for
the future: Archaeology, plants and First Nations’ land use and rights. Pg 86–106 in N. J.
Turner, (Eds.). Plants, people, and places: the roles of ethnobotany and ethnoecology in
Indigenous peoples' land rights in Canada and beyond (Vol. 96). McGill-Queen's PressMQUP.
Lilles, E., Dhar, A., Coates, K. D., & Haeussler, S. (2018). Retention level affects dynamics of
understory plant community recovery in northern temperate hemlock-cedar forests.
Forest Ecology and Management, 421, 3–15.
Lobo, J. M., Jiménez-Valverde, A., & Real, R. (2008). AUC: A misleading measure of the
performance of predictive distribution models. Global Ecology and Biogeography, 17(2),
145–151. https://doi.org/10.1111/j.1466-8238.2007.00358.x
MacKinnon, A., Meidinger, D., & Klinka, K. (1992). Use of the biogeoclimatic ecosystem
classification system in British Columbia. The Forestry Chronicle, 68(1), 100–120.

113

MacKinnon, A., & Pojar, J. (1994). Plants of Coastal British Columbia: Including Washington,
Oregon & Alaska. Lone Pine Publishing.
Mangiafico, S., & Mangiafico, M. S. (2017). Package ‘rcompanion’. 20, 1-71.
Mason, R. (1990). Analytical study of plant/environment interactions in thimbleberry and devil’s
club [MSc Thesis]. University of British Columbia.
Masson-Delmotte, V., P. Zhai, A. Pirani, S..L. Connors, C. Péan, S. Berger, N. Caud, Y. Chen,
L. Goldfarb, M.I. Gomis, M. Huang, K. Leitzell, E. Lonnoy, J.B.R. Matthews, T.K.
Maycock, T. Waterfield, O. Yelekçi, R. Yu, & B. Zhou (eds.)]. (2021). IPCC, 2021:
Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to
the Sixth Assessment Report of the Intergovernmental Panel on Climate Change[
Cambridge, (No. 6). Cambridge University Press, United Kingdom and New York, NY,
USA, In press,. doi:10.1017/9781009157896
McCombs, A. M., & Chittenden, W. W. (1990). The Fraser Valley challenge. Treeline
Publishing,
McCullagh, P., & Nelder, J. A. (1989). Binary data. Generalized linear models.(pp. 98–148)
Springer. USA
McIntire. E. (2001). Within stand properties of natural disturbances in Arrow TSA. End of year
report for Arrow Lakes TSA.
McKinney, M. L., & Lockwood, J. L. (2001). Biotic homogenization: A sequential and selective
process. Biotic homogenization (pp. 1–17). Springer.
McNeill, J. R., & Engelke, P. (2016). The great acceleration: An environmental history of the
Anthropocene since 1945. Harvard University Press.

114

McWethy, D. B., Schoennagel, T., Higuera, P. E., Krawchuk, M., Harvey, B. J., Metcalf, E. C.,
Schultz, C., Miller, C., Metcalf, A. L., & Buma, B. (2019). Rethinking resilience to
wildfire. Nature Sustainability, 2(9), 797–804.
Meidinger, D., & Pojar, J. (1991). Ecosystems of British Columbia. Special Report SeriesMinistry of Forests, British Columbia, 6.
Menzies, C., & Butler, C. (2008). The indigenous foundation of the resource economy of BC’s
north coast. Labour/Le Travailleur, 61, 131–149.
Mishra, S. N., Kumar, D., Kumar, B., & Tiwari, S. (2021). Assessing impact of varying climatic
conditions on distribution of Buchanania cochinchinensis in Jharkhand using species
distribution modeling approach. Current Research in Environmental Sustainability, 3,
100025.
Moola, F. M., & Mallik, A. U. (1998). Morphological plasticity and regeneration strategies of
velvet leaf blueberry (Vaccinium myrtilloides Michx.) following canopy disturbance in
boreal mixedwood forests. Forest Ecology and Management, 111(1), 35–50.
Mucioki, M., Sowerwine, J., Sarna-Wojcicki, D., McCovey, K., & Bourque, S. D. (2022).
Understanding the conservation challenges and needs of culturally significant plant
species through Indigenous Knowledge and species distribution models. Journal for
Nature Conservation, 126285.
Nimasow, G., Nimasow, O. D., Rawat, J. S., Tsering, G., & Litin, T. (2016). Remote sensing and
GIS-based suitability modeling of medicinal plant (Taxus baccata Linn.) in Tawang
district, Arunachal Pradesh, India. Current Science, 219–227.
North, M. E. A., & Teversham, J. M. (1984). The vegetation of the floodplains of the Lower
Fraser, Serpentine and Nicomekl Rivers, 1859 to 1890. Syesis, 17, 47–66.

115

Olden, J. D., Kennard, M. J., Leprieur, F., Tedesco, P. A., Winemiller, K. O., & García-Berthou,
E. (2010). Conservation biogeography of freshwater fishes: Recent progress and future
challenges. Diversity and Distributions, 16(3), 496–513.
Peterson, A. T., & Soberón, J. (2012). Species distribution modeling and ecological niche
modeling: Getting the concepts right. Natureza & Conservação, 10(2), 102–107.
Peterson, A. T., Soberón, J., Pearson, R. G., Anderson, R. P., Martínez-Meyer, E., Nakamura,
M., & Araújo, M. B. (2011). Ecological niches and geographic distributions in R.P.
Anderson, E. Martínez-Meyer, M. Nakamura, M. B Araújo, A. T. Peterson, J. Soberón,
& R. G. Pearson (Eds). Ecological Niches and Geographic Distributions (MPB-49).
Princeton University Press.
Phillips, S. J., Anderson, R. P., & Schapire, R. E. (2006). Maximum entropy modeling of species
geographic distributions. Ecological Modelling, 190(3–4), 231–259.
Pickett, S. T. A., & Kempf, J. S. (1980). Branching patterns in forest shrubs and understory trees
in relation to habitat. New Phytologist, 86(2), 219–228.

Posey, D. A. (1984). A preliminary report on diversified management of tropical forest by the
Kayapo Indians of the Brazilian Amazon. Advances in Economic Botany, 1, 112–126.
Power, S. C., Verboom, G. A., Bond, W. J., Packer, K. F., & Cramer, M. D. (2021). The role of
shade in maintaining alternative stable states between open- and closed-canopy
vegetation. Journal of Ecology, 109(11), 3835–3848. https://doi.org/10.1111/13652745.13760

116

Prevéy, J. S., Parker, L. E., & Harrington, C. A. (2020). Projected impacts of climate change on
the range and phenology of three culturally-important shrub species. PLOS ONE, 15(5),
e0232537. https://doi.org/10.1371/journal.pone.0232537
Prevéy, J. S., Parker, L. E., Harrington, C. A., Lamb, C. T., & Proctor, M. F. (2020). Climate
change shifts in habitat suitability and phenology of huckleberry (Vaccinium
membranaceum). Agricultural and Forest Meteorology, 280, 107803.
https://doi.org/10.1016/j.agrformet.2019.107803
Price, D. T., Alfaro, R. I., Brown, K. J., Flannigan, M. D., Fleming, R. A., Hogg, E. H., Girardin,
M. P., Lakusta, T., Johnston, M., & McKenney, D. W. (2013). Anticipating the
consequences of climate change for Canada’s boreal forest ecosystems. Environmental
Reviews, 21(4), 322–365.
https://www.R-project.org/.
Rabnett, K. (2006). Lower Skeena fish passage assessment highway#16, highway# 37S, and CN
Rail. Skeena Fisheries Commission.
Rabnett, K. (2018). Luutkudziiwus territory current conditions assessment—Interim draft report
Hazelton BC.
Radić, V., Cannon, A. J., Menounos, B., & Gi, N. (2015). Future changes in autumn atmospheric
river events in British Columbia, Canada, as projected by CMIP5 global climate models.
Journal of Geophysical Research: Atmospheres, 120(18), 9279–9302.
Forest Branch. (1930). Report of the Forest Branch. Department of Lands, Victoria, British
Columbia..
Riahi, K., Van Vuuren, D. P., Kriegler, E., Edmonds, J., O’neill, B. C., Fujimori, S., Bauer, N.,
Calvin, K., Dellink, R., & Fricko, O. (2017). The shared socioeconomic pathways and

117

their energy, land use, and greenhouse gas emissions implications: An overview. Global
Environmental Change, 42, 153–168.
Roorbach, A. H. (1999). The ecology of Devil’s club (Oplopanax horridus (JE Smith) Miq.) in
Western Oregon 1999 [MSc Thesis] Oregon State University. Corvallis USA .
Rose, N.A., & Burton, P. J. (2009). Using bioclimatic envelopes to identify temporal corridors in
support of conservation planning in a changing climate. Forest Ecology and
Management, 258, S64–S74.
Rosenvald, R., & Lõhmus, A. (2008). For what, when, and where is green-tree retention better
than clear-cutting? A review of the biodiversity aspects. Forest Ecology and
Management, 255(1), 1–15.
Santini, L., Benítez-López, A., Maiorano, L., Čengić, M., & Huijbregts, M. A. (2021). Assessing
the reliability of species distribution projections in climate change research. Diversity and
Distributions, 27(6), 1035–1050.
Savicky, P. (2015). Spearman’s rank correlation test. Version 0.3–1. Available online:
https://cran. r-project. org/web/packages/pspearman/pspearman. pdf (accessed on 2 May
2022).
Schwab, J. W. (1983). Mass Wasting: October-November 1978 Storm, Rennell Sound, Queen
Charlotte Islands, British Columbia. Research Note-Province of British Columbia,
Ministry of Forests Victoria BC.
Sharma, M., Thakur, R., Sharma, M., Sharma, A. K., & Sharma, A. K. (2020). Changing
Scenario of Medicinal Plants Diversity in Relation to Climate Changes: A Review. Plant
Archives, 20(2), 4389–4400.

118

Siddiqui, M. Z., Chowdhury, A. R., Prasad, N., & Thomas, M. (2014). Buchanania lanzan: A
species of enormous potentials. World Journal of Pharmaceutical Sciences, 374–379.
Siegel, A. F. (1982). Robust regression using repeated medians. Biometrika, 69(1), 242–244.
Sinclair, S. J., White, M. D., & Newell, G. R. (2010). How useful are species distribution models
for managing biodiversity under future climates? Ecology and Society, 15(1).
Https://www.jstor.org/stable/26268111
Smith, H. I.(1927a).: The Carrier Indians of British Columbia. National Museum of Canada.
Ottawa.
Smith, H. I.. (1927b). The Bella Coola Indians of British Columbia. National Museum of
Canada. Ottawa
Smith, H. I. (1927c). The Tsimshian (Gitksan) Indians of British Columbia. National Museum of
Canada. Ottawa
Swart, N. C., Cole, J. N., Kharin, V. V., Lazare, M., Scinocca, J. F., Gillett, N. P., ... & Winter,
B. (2019). The Canadian earth system model version 5 (CanESM5). Geoscientific Model
Development, 12(11), 4823-4873.
Tai, J., Cheung, S. S. C., Ou, D., Warnock, G. L., & Hasman, D. (2014). Antiproliferation
activity of Devil’s club (Oplopanax horridus) and anticancer agents on human pancreatic
cancer multicellular spheroids. Phytomedicine, 21(4), 506–514.
Thompson, V., Kennedy-Asser, A. T., Vosper, E., Lo, Y. E., Huntingford, C., Andrews, O.,
Collins, M., Hegerl, G. C., & Mitchell, D. (2022). The 2021 western North America heat
wave among the most extreme events ever recorded globally. Science Advances, 8(18),
eabm6860.

119

Thuiller, W. (2004). Patterns and uncertainties of species’ range shifts under climate change.
Global Change Biology, 10(12), 2020–2027.
Thuiller, W., Lavorel, S., & Araújo, M. B. (2005). Niche properties and geographical extent as
predictors of species sensitivity to climate change. Global Ecology and Biogeography,
14(4), 347–357.
Tilford, G. L. (1997). Edible and medicinal plants of the west. Mountain Press Publishing.
Turner, N. (1982).. Traditional use of devil’s-club (Oplopanax horridus; Araliaceae) by native
peoples in western North America. Journal of Ethnobiology, 2(1), 17-38.
Turner, N. (2014). Ancient pathways, ancestral knowledge: Ethnobotany and ecological wisdom
of indigenous peoples of northwestern North America. McGill-Queen’s Press
Turner, N. J., & Berkes, F. (2006). Coming to understanding: Developing conservation through
incremental learning in the Pacific Northwest. Human Ecology, 34(4), 495–513.
Turner, N. J., Davidson, F. E., & Enrico, J. J. (2004). Plants of Haida Gwaii. Winlaw, BC: Sono
Nis Press.
Turner, N. J., Deur, D., & Mellott, C. R. (2011). “Up on the mountain”: Ethnobotanical
importance of montane sites In Pacific coastal North America. Journal of Ethnobiology,
31(1), 4–43.
United States Geological Survey (2020) (USGS) Landsat 8 Image Archive WRS path 52 Row 21
Uprety, Y., Asselin, H., Dhakal, A., & Julien, N. (2012). Traditional use of medicinal plants in
the boreal forest of Canada: Review and perspectives. Journal of Ethnobiology and
Ethnomedicine, 8(1), 7. https://doi.org/10.1186/1746-4269-8-7

120

Valladares, F., & Niinemets, Ü. (2008). Shade Tolerance, a Key Plant Feature of Complex
Nature and Consequences. Annual Review of Ecology, Evolution, and Systematics, 39,
237–257.
Vapnik, V. (1998). The support vector method of function estimation. Nonlinear modeling:
Advanced black-box techniques (pp. 55-85). Boston, MA: Springer USA.
Vuksan, V., Xu, Z. Z., Jovanovski, E., Jenkins, A. L., Beljan-Zdravkovic, U., Sievenpiper, J. L.,
Mark Stavro, P., Zurbau, A., Duvnjak, L., & Li, M. Z. (2019). Efficacy and safety of
American ginseng (Panax quinquefolius L.) extract on glycemic control and
cardiovascular risk factors in individuals with type 2 diabetes: A double-blind,
randomized, cross-over clinical trial. European Journal of Nutrition, 58(3), 1237–1245.
Wan, J., Wang, C., Yu, J., Nie, S., Han, S., Liu, J., Zu, Y., & Wang, Q. (2016). Developing
conservation strategies for Pinus koraiensis and Eleutherococcus senticosus by using
model-based geographic distributions. Journal of Forestry Research, 27(2), 389–400.
Wang, C.Z., Zhang, Z., Huang, W.-H., Du, G.J., Wen, X.-D., Calway, T., Yu, C., Nass, R., Zhao,
J., & Du, W. (2013). Identification of potential anticancer compounds from Oplopanax
horridus. Phytomedicine, 20(11), 999–1006.
Wang, J., Shao, L., Wang, C.Z., Zhou, H.-H., Yuan, C.-S., & Huang, W.H. (2019). Synergetic
inhibition of human colorectal cancer cells by combining polyyne-enriched fraction from
Oplopanax elatus and irinotecan. Nutrition and Cancer, 71(3), 472–482.
Warman, R. D. (2014). Global wood production from natural forests has peaked. Biodiversity
and Conservation, 23(5), 1063–1078.
Wedley, J. R. (1990). Laying the Golden Egg: The Coalition Government’s Role in Post-War
Northern Development. BC Studies: The British Columbian Quarterly, 88, 58–92.

121

Wong, C., Dorner, B., & Sandmann, H. (2003). Estimating historical variability of natural
disturbances in British Columbia. Land Management Handbook. Ministry of Forests
Forest Science Program Ministry of Sustainable Resource Management Resource
Planning Branch. Victoria BC.
Wong, I. W., Bloom, R., McNicol, D. K., Fong, P., Russell, R., & Chen, X. (2007). Species at
risk: Data and knowledge management within the WILDSPACETM Decision Support
System. Environmental Modelling & Software, 22(4), 423–430.
WorldClim (2021) available at https://www.worldclim.org .accessed on January 1st 2021
Yang, J., Lee, Y.J., Kwon, Y. S., & Kim, M. J. (2018). Anticancer activity of an Oplopanax
elatus stem extract and biologically active isolated compounds. Current Pharmaceutical
Biotechnology, 19(3), 258–264.
Zhang, Q., Wen, J., Chang, Z. Q., Xie, C. X., & Song, J. Y. (2018). Evaluation and prediction of
ecological suitability of medicinal plant American ginseng (Panax
quinquefolius). Chinese Herbal Medicines, 10(1), 80-85.
Zhang, X., & World Health Organization, (2002). Traditional medicine strategy 2002 2005.

122

Appendix
Appendix 1.
List of Sampling Site Locations

1

1970

1970-01

544

1970

55°17'48"N

127°16'21"W

2

1970

1970-02

701

1976

55°18'0"N

127°17'29"W

3

1970

1970-03

616

1978

55°19'0"N

127°16'0" W

4

1970

1970-04

610

1978

55°19'38"N

127°14'46"W

5

1970

1970-05

988

1979

55°26'0"N

127°10'23"W

6

1980

1980-01

672

1984

55°17'23"N

127°18'0"W

7

1980

1980-02

766

1980

55°19'53"N

127°13'22"W

8

1980

1980-03

760

1983

55°19'49"N

127°12'1"W

9

1980

1980-04

895

1986

55°19'52"N

127°9'32"W

10

1990

1990-01

843

1996

55°20'30"N

127°20'47"W

11

1990

1990-02

816

1993

55°20'46"N

127°20'5"W

12

1990

1990-03

908

1994

55°21'4"N

127°22'11"W

13

1990

1990-04

902

1998

55°21'25"N

127°22'45"W

14

2000

2000-01

668

2012

55°17'51"N

127°17'22"W

15

2000

2000-02

941

2000

55°21'19"N

127°20'14"W

16

2000

2000-03

968

2012

55°25'50"N

127°9'57"W

17

control

3000-01

562

1739

55°18'52"N

127°16'35"W

18

control

3000-02

686

1705

55°18'59"N

127°17'30"W

19

control

3000-03

726

1747

55°19'57"N

127°19'14"W

20

control

3000-04

742

1884

55°19'58"N

127°12'49"W

21

trad-use

4000-01

474

1847

55°17'18"N

127°11'20"W

22

trad-use

4000-02

401

1892

55°16'34"N

127°18'4"W

23

trad-use

4000-03

411

1892

55°16'38"N

127°17'58"W

24

trad-use

4000-04

515

1887

55°17'35"N

127°12'37"W

123

Appendix 2.
Explanatory list of BIOCLIM ecological variables Bioclimatic variables.
Bioclimatic variables are derived from the monthly temperature and rainfall values in order to
generate more biologically meaningful variables. These are often used in species distribution
modeling and related ecological modeling techniques. The bioclimatic variables represent annual
trends (e.g., mean annual temperature, annual precipitation) seasonality (e.g., annual range in
temperature and precipitation) and extreme or limiting environmental factors (e.g., temperature
of the coldest and warmest month, and precipitation of the wet and dry quarters). A quarter is a
period of three months (quarter of the year). (WorldClim.org 2022)
They are coded as follows:
BIO1
BIO2
BIO3
BIO4
BIO5
BIO6
BIO7
BIO8
BIO9
BIO10
BIO11
BIO12
BIO13
BIO14
BIO15
BIO16
BIO17
BIO18
BIO19

Annual Mean Temperature
Mean Diurnal Range (Mean of monthly (max temp - min temp)
Isothermality (BIO2/BIO7) (×100)
Temperature Seasonality (standard deviation ×100)
Max Temperature of Warmest Month
Min Temperature of Coldest Month
Temperature Annual Range (BIO5-BIO6)
Mean Temperature of Wettest Quarter
Mean Temperature of Driest Quarter
Mean Temperature of Warmest Quarter
Mean Temperature of Coldest Quarter
Annual Precipitation
Precipitation of Wettest Month
Precipitation of Driest Month
Precipitation Seasonality (Coefficient of Variation)
Precipitation of Wettest Quarter
Precipitation of Driest Quarter
Precipitation of Warmest Quarter
Precipitation of Coldest Quarter

124

Appendix 3.
Table showing the results of the Shapiro-Wilk normality test from the values provided from the
field sampling data.

Variable
Stem Length
Stem Abundance
Stem Diameter
Elevation
Shade
Leaf Area

Shapiro-Wilk normality test
W - Value
0.95881
0.8881
0.97787
0.97519
0.93147
0.19857

p-value
0.01776
1.55E-06
0.2273
0.08917
0.0001693
< 2.2e-16

125

Appendix 4.
Model Evaluation based on AUC scores for Domain, GLM, Gaussian, MaxEnt and SVM model
output.
DOMAIN
Class
n presences
n absences
AUC
cor
max TPR+TNR at

Model Evaluation
18
144
0.8572531
0.4198339
0.2583034

GENERAL LINEAR REGRESSION
Class
Model Evaluation
n presences
18
n absences
144
AUC
0.8996914
cor
0.3854456
max TPR+TNR at
-1.896563
GAUSSIAN
class
n presences
n absences
AUC
cor
max TPR+TNR at

Model Evaluation
18
144
0.8695988
0.4407287
0.1809112

SUPPORT VECTOR MACHINES
class
Model Evaluation
n presences
18
n absences
144
AUC
0.9205247
cor
0.5239749
max TPR+TNR at
0.03078363
MAXENT
class
n presences
n absences
AUC
cor
max TPR+TNR at

Model Evaluation
18
144
0.9085648
0.546431
0.4980119

126

Appendix 5.
Weighted Mean of Suitability map

Map showing areas weighted mean of GLM, BIOCLIM and Support Vector Machine models.
The greener colour indicates the most suitable areas and grey area means least suitable area for
devil’s club habitat.

127

Appendix 6
Ecological variable contribution plot for predicted habitat in MaxEnt model for British Columbia

128

Appendix 7
Maximum Entropy Response plots for devil’s club for British Columbia predictive habitat
models

129