Plant succession on glacial moraines in the Arctic Brooks Range along a >125,000-year glacial chronosequence/toposequence

ABSTRACT

Widespread glacial retreat is now occurring in many arctic mountain ranges, yet little is known about primary succession following deglaciation in these settings. Newly created habitats could provide refugia for flora and fauna whose ranges are threatened elsewhere by rapid warming. To assess vegetation responses to glacial retreat in an arctic–alpine setting, we first describe plant community development on two recently deglaciated moraines in the Brooks Range. We then compare these recent communities with communities developed along a moraine chronosequence that spans >125,000 years and ranges in altitude between 800 and 1,700 m.a.s.l. Results show that (1) within twenty-two to thirty-six years following deglaciation, primary succession begins with the assembly of small communities of eight to thirteen vascular and nonvascular plant species; (2) species turnover is low, with many pioneer taxa, particularly lichens, persisting at the oldest sites and across all altitudes; and (3) overall, succession is directional and slow, with species richness increasing for up to 25,000 years, and percentage vegetation cover reaching >100 percent on the oldest glacial deposits. This is the first vegetation study on primary succession in the high central Brooks Range, and it supplies a previously missing alpine element within a vegetation transect across northern Alaska’s bioclimatic gradient.

KEYWORDS:

• Lichenometry
• primary succession
• deglaciation
• Arctic
• alpine

Introduction

Globally, glaciers have been retreating since the end of the Little Ice Age (ca. 1850 CE), in large part due to climate change (Schumann, Gewolf, and Tackenberg 2016). Glaciers remaining in the highest peaks of the Brooks Range have been retreating rapidly over the past several centuries and are now restricted to high north-facing cirques (Ellis 1978; Calkin and Ellis 1980). As these small glaciers disappear and new surfaces are exposed, fresh substrates become available for the establishment of pioneer plant communities. Pioneer plant communities are a key component of plant community development in arctic environments (Mori et al. 2008; Monzón, Moyer-Horner, and Palamar 2011). Primary succession affects soil development, nutrient cycling, and plant community assembly (Fastie 1995). In the rapidly changing Arctic, understanding the trajectories and rates of these ecological shifts can help predict and manage the ecological impacts of ongoing climate changes and other disturbances (Post et al. 2009). Pioneer plant communities also create new habitats that may serve as refugia for flora and fauna whose ranges are being shifted by polar amplification of climate change (Chang and Turner 2019).

Previous studies have described plant communities at lower altitudes in the northern foothills of the Brooks Range (Jorgenson 1984; M. D. Walker, Walker, and Auerbach 1994; D. A. Walker and Walker 1996; D. A. Walker et al. 2014). Pollen records have been used to interpret ecological legacies of Itkillik (late Pleistocene) and Sagavanirktok (middle Pleistocene) glaciations (Oswald et al. 2003). However, little attention has been paid to ecological succession associated with recent deglaciation at high elevations in the Alaskan Arctic. Rapid glacial retreat in the central Brooks Range creates an ideal setting to study the patterns and processes of primary succession in an arctic–alpine environment and how they relate to older, established plant communities on similar glacial deposits. Linking the structure and development of pioneer and early-successional communities to that of late-successional and mature communities is a missing component necessary for managing and predicting changes occurring in fragile and globally important arctic ecosystems.

Processes of primary succession have been extensively documented in temperate, boreal, subarctic, and alpine environments (Cox 1933; Matthews and Whittaker 1987; Chapin et al. 1994; Fastie 1995; C. C. Jones and Del Moral 2005; Hollingsworth et al. 2010; Carlson et al. 2014; Schumann, Gewolf, and Tackenberg 2016; Matthews et al. 2018) and some arctic environments (Svoboda and Henry 1987; Bliss and Gold 1994; G. A. Jones and Henry 2003; Mori et al. 2008; Mori, Uchida, and Kanda 2013), but relatively little is known about the earliest stages of plant succession that immediately follow deglaciation in arctic mountain ranges, in part due to the remoteness of these sites. Perhaps the most thorough recent study is that of two contrasting high Arctic glacial foreland chronosequences up to 2,000 years old in Spitsbergen (Hodkinson, Coulson, and Webb 2003). However, these moraines are at sea level and are strongly influenced by the maritime climate, as are most previously studied arctic–alpine, primary succession study sites (Bliss and Gold 1994; Hodkinson, Coulson, and Webb 2003; G. A. Jones and Henry 2003). In Spitsbergen, cyanobacteria were dominant pioneer colonizers of soil surfaces for the first 60 years of primary succession, vascular plants increased in species richness for the first 100 years, and shrubs slowly became dominant with time (Hodkinson, Coulson, and Webb 2003). Development of soil and plant communities occurred across different timelines, depending on the availability of nutrients and variability of microclimates. Recruitment of pioneer species was limited to specific pools of species able to disperse to freshly exposed substrates and amenable to the specific ecological conditions available on each site (Hodkinson et al. 2003).

Although there are many commonalities between the Svalbard glacial forelands and forelands in the Brooks Range, such as the arctic climate and limitations in species dispersal and recruitment, there are important environmental differences between sites that likely affect pioneer community establishment and subsequent plant community development. There are major environmental differences between coastal glacial forelands, like those in Svalbard, and the more continental alpine setting of the Grizzly Glacier. Based on previously described factors of soil quality, pH, moisture, and climate as being important driving factors of pioneer plant community establishment and successional processes (Svoboda and Henry 1987; Chapin et al. 1994; Hodkinson, Coulson, and Webb 2003; Robbins and Matthews 2009), primary succession in interior arctic locations generally occur at a slower pace and with a different suite of pioneer species, relative to other arctic–alpine successional processes (Svoboda and Henry 1987; Bliss and Gold 1994; Hodkinson, Coulson, and Webb 2003; G. A. Jones and Henry 2003).

This study has two main objectives: (1) establish the timeline of moraine stabilization and plant succession from the high north-facing Grizzly Glacier cirque, downslope across the northern foothills of the Brooks Range, and (2) document how plant communities change and persist relative to environmental factors that vary along gradients of time since deglaciation and altitude and identify the factors most important in directing plant community structure and development over the course of succession. Specific research questions include (1) How quickly does succession begin following deglaciation, and how long does it continue afterwards? (2) What is the species composition and structure of cirque pioneer communities, and what are the similarities and differences relative to other successional arctic moraine communities? (3) Does succession progress directionally in this system, and how much turnover occurs in plant species and growth forms during succession?

Materials and methods

Study area

The Brooks Range extends 1,000 km from the Chukchi Sea to the Alaska/Yukon border and divides boreal forest to its south from tundra to its north. The Brooks Range has been glaciated repeatedly over the last 3 million years (Hamilton and Porter 1975; D. A. Walker et al. 2014; Pendleton et al. 2015). Here we utilize a chronosequence consisting of six glacial deposits whose minimum ages range from 22 years to approximately 125,000 years since deglaciation (Table 1; Hamilton 2003). These moraines are located along an altitudinal gradient extending from the crest of the Brooks Range near Atigun Pass, northward into the Southern Arctic Foothills (Figure 1). Study sites on these glacial deposits were accessed on foot and occasionally by helicopter from the Dalton Highway and the Toolik Research Station. We studied six glacial deposits: (1) the End of the Little Ice Age moraine in the Grizzly cirque (ELIA, 22–36 years since deglaciation), (2) the neoglacial moraine in the Grizzly cirque (NEO, 238–500 years since deglaciation; Ellis 1978), (3) an Itkillik II alpine glacial deposit near Atigun Pass (ITKII-AP, 10,000–25,000; Hamilton and Porter 1975; Pendleton et al. 2015; Jim Ellis, Ellis Geospatial, Cambria, NY, personal communication, 2017), (4) an Itkillik II foothill moraine (ITKII, 10,000–25,000 years since deglaciation; Hamilton 2003), (5) an Itkillik I foothill moraine (ITKI, 25,000–75,000 years since deglaciation; Hamilton 2003), and (6) a Sagavanirktok foothill moraine (SAG, 125,000–780,000 years since deglaciation; Hamilton 2003). Three of these glacial deposits are located in the alpine above 900 m elevation (ELIA, NEO, and ITKII-AP) near Atigun Pass, and three are in the Arctic foothills northeast of Toolik Lake (ITKII-N, ITKI, SAG; Figure 1, Table 1). The two oldest glacial deposits only occur in the foothills, and the two youngest glacial deposits only occur near the range crest. The Itkillik II (ITKII) glacial deposits occur in both the foothills and the alpine (Hamilton and Porter 1975; Pendleton et al. 2015; Jim Ellis, Ellis Geospatial, Cambria, NY, personal communication, 2017) and are used here to also assess effects of elevation on plant community development on glacial deposits near the same age—one near the summit of Atigun Pass (ITKII-AP; ~1,500 m.a.s.l.) and one in the foothills near Slope Mountain (ITKII-N; ~800 m.a.s.l.; Figure 1, Table 1).

Figure 1. Study site locations. From youngest to oldest, the Grizzly Glacier site includes ELIA and NEO, Itkillik II aged glacial deposits are located at ITKII-AP and near the toe of Slope Mountain in the foothills (ITKII-N), and ITKI and SAG are located between Imnavait Creek and Toolik Lake in the foothills. Map created by Dr. Martha Raynolds at the University of Alaska Fairbanks.

Table 1. Study site characteristics and plots sampled per site.

Glacial deposit	Average elevation (M)	Latitude and longitude (DD)	Estimated age (years BP)	Rock relevés (N)	Fine-grained relevés (N)
ELIA	1,700	68.1133 N, 149.4426 W	22–36	3	3
NEO	1,700	68.1132 N, 149.4386 W	238–500	3	3
ITKII-AP	1,500	68.1299 N, 149.4770 W	10,000–15,000	4	4
ITKII-N	800	68.6762 N, 149.0888 W	10,000–15,000	5	5
ITKI	900	68.8672 N, 149.5101 W	25,000–75,000	3	3
SAG	950	68.6125 N, 149.30016 W	125,000–150,000	3	3
(N plots = 42)	Δaltitude = 900 m		$\mathrm{\Delta }\mathit{time}=\ge$124,978 years	N = 21	N = 21

In the late 1970s, Parker Calkin and Jim Ellis conducted a series of glacial surveys across central Brooks Range cirques (Figure 2). The colonization of recently deposited moraine boulders by common pioneer lichen species was used to create a lichenometric age curve for estimating the ages of moraines (Ellis 1978; Calkin and Ellis 1980, 1981). These lichenometric surveys inspired this study and provide the basis for comparing lichenometric measurements of more recently deposited moraines. We used Calkin and Ellis’s methodology to estimate the time since deglaciation of the most recently deposited moraine in the Grizzly glacial cirque, henceforth referred to as the End of the Little Ice Age moraine (Figure 2), and the likely successional order of pioneer species colonizing both the neoglacial and ELIA moraines in this cirque.

Figure 2. (Top) Grizzly Glacier cirque photographed and used with permission by Jim Ellis in July 1977. (Bottom) The same general scene photographed by SK (formerly S. Gowan) in July 2017. These photographs were taken from slightly different locations and distances from the glacier. Red arrows indicate moraine features for comparison of ice mass.

There have been few studies of the vegetation in the alpine zone of the central Brooks Range, despite the presence of the Dalton Highway, which provides access to the peaks in the vicinity of the Atigun Pass. A brief reconnaissance survey and map of the vegetation in the immediate vicinity of the summit of Atigun pass was made during 1975–1976 as part of a survey of vegetation along the Dalton Highway prepared for the U.S. Army Cold Regions Research and Engineering Laboratory (P. Webber et al. 1979). The most thorough vegetation study of any alpine area in the Brooks Range was that of David Cooper (Cooper 1986, 1989a, 1989b), who studied the vegetation of the Arrigetch Peaks, approximately 210 km west of the Grizzly Glacier on the south side of the Brooks Range. Cooper sampled 439 plots in forty-nine different plant associations, including several on glacial moraines.

The Grizzly Glacier, would fall into bioclimatic elevation belt A based on its elevation of approximately 1,700 m within Bioclimate Subzone E of the Circumpolar Arctic Vegetation Map (CAVM Mapping Team et al. 2003). However, the summer warmth index (sum of the monthly mean temperatures above 0°C) at the nearby Atigun Pass Snotel site (elevation 1,453 m) averaged 17.8°C during 2008–2018, which would place Atigun Pass in elevation belt D based on current summer temperatures and the base of Grizzly Glacier, which is approximately 250 m higher, marginally in elevation belt C. The lower-elevation foothill deposits broadly fit into bioclimatic subzone E, with June, July, August mean monthly temperatures of 9°C to 13°C and summer warmth index averaging 33.1°C (CAVM Mapping Team et al. 2003).

Site selection and field sampling

We measured a suite of ecological and environmental parameters known to impact plant community structure within forty-two 2 m × 2 m plots. Morainal ridges from the high alpine to the foothills share a similar geomorphology; namely, a mixture of boulder-covered areas and areas of fine-grained soil, which differ substantially in the type of plant communities they support. Thus, it was important to survey plant communities on both substrates to obtain a picture of plant succession on moraine ridges. Three to five plots were placed on both “rocky” deposits and “fine-grained” sediment at least 10 m apart on each glacial surface (Figure 3). Plot placement was subjective and used “centralized replicate” sampling of the Braun-Blanquet vegetation classification approach (Westhoff and van der Maarel 1978). We selected relatively large homogeneous flat or gently sloping “rocky” sites, dominated by well-drained, large boulders with minimal bare soil (Figure 3a), and corresponding “fine-grain” sites on the crests of glacial deposits where boulders were scarce (Figure 3b). Once sites with rocky and fine-grained characteristics were identified on each glacial deposit, plots were placed centrally within large areas of both substrates to avoid edge effects. Large homogeneous areas were less common on the youngest site, and three plots of each substrate type was the maximum number possible on the youngest (ELIA and NEO) moraines in the Grizzly glacier cirque. We maintained a replicate number of at least three plots per substrate consistent on ELIA, NEO, ITKI, and SAG glacial surfaces, but on ITKII aged deposits, which occurred at both high and low elevations, five plots were placed to more thoroughly compare plant communities on glacial deposits of roughly the same age at different elevations.

Figure 3. Examples of (A) rock and (B) fine-grained substrate plots on each glacial deposit surface in order of the estimated age since deglaciation.

Plant communities were sampled in the field according to the Braun-Blanquet approach (Westhoff and Maarel 1978) as modified for arctic plant communities (D. A. Walker et al. 2016). Cover of each vascular plant, lichen, and bryophyte was visually estimated using Braun-Blanquet cover abundance scores (r = rare, + = common but <1 percent, 1 = 1–5 percent, 2 = 6–25 percent, 3 = 26–50 percent, 4 = 51–75 percent, 5 = 76–100 percent; Barkman, Moravec, and Rauschert 1976) and numerically transformed for use in PC-ORD software (McCune and Mefford 2011; Supplementary Table A). Plots were established at least 1 km from the Dalton Highway when possible, on distinct glacial deposit ridges previously identified and dated, with the exception of the youngest moraine (ELIA) dated in this study (Calkin and Ellis 1980; Hamilton 2003). The ITKII-AP alpine site is less than 1 km from the road in places, and this site has not been mapped as a distinct glacial moraine; however, based on descriptions in the literature, the ITKII-AP feature was likely deposited during the Itkillik II glaciation phase (10,000–25,000 years ago; Hamilton and Porter 1975; Hamilton 2003; Pendleton et al. 2015; Jim Ellis, Ellis Geospatial, Cambria, NY, personal communication, 2017). All sampling sites were established at locations where siliceous rather than carbonate rock types predominate.

Plant collection, curation, and identification

Vascular plant collections were identified by the authors with aid from Carolyn Parker, taxonomic research professional at the University of Alaska Fairbanks Herbarium. Nomenclature for vascular plants was originally sourced from Flora of Alaska and Neighboring Territories (Hultén 1990) and for bryophytes from Flora of North America, Bryophyta (Miller and Miller 2007). Nomenclature was then aligned with the latest available version of Flora of North America respective volumes using the online database (http://floranorthamerica.org) and were synchronized with the Panarctic flora checklist when possible (http://panarcticflora.org/). Nomenclature for lichens follows the Consortium of North American Lichen Herbaria (lichenportal.org) and American Arctic Lichens Volumes 1 and 2 (Thomson and Brehmer 1984). Lichen species were identified by S.K. using standard macroscopic field identification, chemical spot tests, and microscopy. Lichen identification was verified by Dr. Mikhail P. Zhurbenko, lead researcher of the Systematics and Geography of Fungi Laboratory at the Komarov Botanical Institute Russian Academy of Sciences. Voucher specimens were curated in the Komarov Botanical Institute nonvascular herbarium. Bryophyte collections were verified and curated by David Toren of the California Academy of Sciences Herbarium. A subset of bryophyte samples was verified by Andrew Cortese at the University of Alaska Fairbanks using standard microscopy techniques.

A representative sample of every plant species found on each plot was collected for verification and comparison purposes whenever possible. Vascular plants were stored in a plant press while in the field, and cryptogams were stored either together by plot or separately in labeled wax, paper, or plastic bags and allowed to air dry. In the laboratory, plants were sorted, identified, and given a collection number for further verification purposes.

Environmental data

To compare seasonal changes in temperature at each site over the course of one year, a temperature logging device (iButton; https://www.maximintegrated.com) was placed in the southwest corner of each plot just below the litter/moss layer of the soil from August 2017 through July 2018. If the litter/moss layer was absent in a given plot (e.g., on rocky substrate), the iButton was placed in a shallow depression and covered with a stone typical of those present in the plot. Using a combination of data gathered via iButtons and data from two weather stations located at Imnavait creek (SNOTEL #968) and Atigun Pass (SNOTEL # 957), summer warmth index (the sum of mean monthly temperatures >0°C), mean annual freeze season (months with an average temperature below 0°C), and mean annual thaw season (months with an average temperature above 0°C) temperatures were calculated for the years 2006 to 2018.

Where possible, 235 cm³ of soil was collected from below the litter layer to a depth of approximately 10 cm from immediately outside of each plot. In the laboratory, moisture content was measured first by comparing wet and dry weight of each sample for gravimetric soil moisture and volumetric soil moisture when it was possible to collect an intact sample with the soil can (Peters 1965). Wet and dry soil color and chroma were determined using a Munsell color book (Munsell Color Firm 2009). Soil particle size was analyzed by sieving the soil through a 2-mm screen to remove gravel. Percentage gravel was determined by subtracting the weight of the gravel removed via the sieve from the total dry weight of the soil sample. The hydrometer method was used to determine percent sand, silt, and clay (Bouyoucos 1962). pH was measured with a glass electrode pH meter and the paste method (McLean 1982). Additionally, eighteen site characteristics (slope, aspect, elevation, landform, surficial geology, surficial geomorphology, microsite, site moisture, soil moisture, soil texture of top mineral horizon, glacial geology, topographic position, habitat type, estimated snow duration, disturbance degree, disturbance type, physical stability, and exposure) together with the percentage cover of all plant growth forms, soil, rock, and water were estimated and recorded for each plot in the field (see Supplementary Table B for the complete environmental matrix).

Data analysis

Two matrices of species cover abundance and corresponding environmental data were compiled by plot for use in the ordination program PC-ORD 6 (McCune and Mefford 2011). The species matrix contains plot numbers (columns) and species names (rows) with Braun-Blanquet cover abundance values in each cell (Supplementary Table A). The environmental matrix contains plot numbers (columns) and environmental information (rows; Supplementary Table B). Vegetation and environmental data were deposited in the Arctic Vegetation Archive database at the Alaska Geobotany Center located on the University of Alaska Fairbanks campus (https://www.geobotany.uaf.edu/ava/) and in the Arctic Data Center database (Kasanke et al. 2023).

Lichenometry

Lichenometric measurements were used to estimate the ages of rock surfaces exposed since deglaciation based on the diameters of lichen thalli whose growth rates are known (Andrews and Webber 1969; Loso and Doak 2006; Wiles, Barclay, and Young 2010; Beschel 1961). In the late 1970s, Calkin and Ellis (1980) constructed a lichen growth curve for the Central Brooks Range region, which they used to date the stabilization of the NEO moraine in the Grizzly Glacier cirque at approximately 500 years old. Calkin and Ellis established a baseline growth curve for three lichen taxa by measuring the diameter of the largest lichen thalli found on mine tailings and gravestones of known age in the region, bolstering these measurements with radiocarbon and dendrochronological data when possible (Calkin and Ellis 1980, 1981). When a glacier retreats and exposes a moraine, and when this moraine is no longer undergoing thermokarsting caused by the melting of buried glacial ice, large rocks stabilize and allow the colonization of crustose lichens; this is when a moraine is considered to be stabilized for the purpose of this study. Following the methods of Calkin and Ellis (1980), lichenometry was used to estimate stabilization of the ELIA moraine and general rates of pioneer lichen colonization and to verify the estimated age since deglaciation of the NEO moraine as dated by Calkin and Ellis.

Rhizocarpon geographicum s.l. is an easily recognized saxicolous crustose lichen commonly used for lichenometry (P. J. Webber and Andrews 1973; Calkin and Ellis 1980; Innes 1985; Wiles, Barclay, and Young 2010). It grows on a variety of rock substrates and environments, including arctic–alpine. It should be noted that here the abbreviation s.l. (sensu lato) is used to designate a complex of similar Rhizocarpon species that share similar growth rates in similar environments (Ellis 1978; Calkin and Ellis 1980; Innes 1985; Wiles, Barclay, and Young 2010). R. geographicum s.l. has been shown to grow at a rate of about 3 mm/100 years in the central Brooks Range following an initial “great growth period” lasting about 200 years when it grows at about 10 mm/100 years, (Calkin and Ellis 1980; Bradwell and Armstrong 2007).

We also used faster growing lichen species in the genera Pseudephebe (P. minuscula and P. pubescens) and Umbilicaria (U. proboscidea, U. hyperborea, and U. cylindrica) to help estimate the time of moraine stabilization in the Grizzly glacier cirque (Andrews and Webber 1969; Calkin and Ellis 1980). These groups of species within the genera Pseudephebe and Umbilicaria share similar growth rates (P. J. Webber and Andrews 1973; Calkin and Ellis 1980) and will be discussed together as Pseudephebe and Umbilicaria spp. henceforth. Pseudephebe and Umbilicaria spp. grow approximately seven times faster than R. geographicum and do not live as long (hundreds vs. thousands of years; Calkin and Ellis 1980); they are also more sensitive to disturbance and thus cannot successfully colonize a surface until it is completely stable (Innes 1985; Nash 1996).

To improve the accuracy of our age estimates of moraine stabilization and to account for the possibility of “inherited” thalli (lichens that were already growing on boulders that fell from the headwall of the cirque during deglaciation), we used a combination of the “five largest thalli” averaging approach (Innes 1985; Haworth, Calkin, and Ellis 1986; McCarthy 2021) and the “single largest thallus” approach (Calkin and Ellis 1980; Haworth, Calkin, and Ellis 1986; Solomina and Calkin 2003). The largest thalli of the target lichen groups were measured on multiple boulders, across both the NEO and ELIA moraines in the Grizzly Glacier cirque. This was accomplished by surveying the east-facing, west-facing, and north-facing lateral moraines on foot with a team of three field technicians who measured the largest thalli of both groups of lichen taxa whenever a colonized boulder was located. A minimum of ten measurements of each targeted lichen group were taken on each lateral moraine. From this data set, the sizes of the five largest lichen thalli were averaged to estimate stabilization dates of each moraine. The five largest R. geographicum s.l. thalli on each moraine were used to calculate for the upper age estimate boundary, and the five largest Pseudephebe spp. and Umbilicaria spp. on each moraine were used to calculate the lower age estimate boundary. Age estimations were calculated according to Calkin and Ellis’s lichenometric growth curve, adjusting for growth rate changes in both groups of lichens over time (Calkin and Ellis 1980). Lichenometric measurements were taken on all glacial surfaces across the sequence, but measurements on moraines >500 years old are not reported here because it has been shown that the accuracy of lichenometry degrades over time due to the effects of rock weathering, thallus aging, and turnover within the lichen population (Innes 1985; Loso and Doak 2006; McCarthy 2021).

Cluster analysis

A hierarchical dendrogram approach was used to assign plots to discrete groups based on similarities in plant community composition (PC-ORD 6; McCune and Mefford 2011). Sørenson’s distance measurement and a flexible beta group linkage method ($\mathit{\beta }=0.25$) were selected as analytical parameters to create the dendrogram. This resulted in seven major clusters of plant community types. Each community type according to the cluster analysis is described in detail in Appendix A according to the Braun-Blanquet approach (Westhoff and Maarel 1978), and community types were defined in accordance with descriptions in Mucina et al.'s Vegetation of Europe (Mucina et al. 2016), which includes arctic and alpine community types.

Synoptic table and indicator species analysis

An indicator species analysis was used to define dominant and characteristic species of each community type identified in the cluster analysis (Chytrý et al. 2002; PC-ORD 6 indicator species analysis tool; McCune and Mefford 2011). Characteristic species of each community type were identified and ranked based on fidelity scores. Fidelity is a statistical measurement of the concentration of species occurrences within vegetation units (Chytrý et al. 2002). The fidelity of a certain plant to a given community type is measured by a phi value. The phi value ranges from zero to one, with values closer to one indicating a higher degree of fidelity (Chytrý et al. 2002). Species with phi coefficients ≥0.4 were considered somewhat characteristic, those with phi coefficients ≥0.5 were considered characteristic and those with phi coefficients ≥0.8 were considered highly characteristic. Because the phi value has no test of statistical significance associated with it (Chytry et al. 2002), a Monte Carlo randomization test was used to assess whether the phi values associated with each species within a cluster was less than 95 percent likely to have been assigned by chance (p >.05; 4,999 randomized runs). A synoptic table was then prepared by hand, combining the results of the cluster analysis and the indicator species analysis (Appendix B). For each cluster in the synoptic table, characteristic species were organized by fidelity (phi values from the indicator species analysis) from lowest to highest. The frequency of each species’ occurrence in each cluster was determined by averaging species frequency across all plots included in each cluster (Appendix B). Taxa that were noncharacteristic according to the indicator species analysis were removed from the clusters and included in the synoptic table section labeled noncharacteristic taxa (Appendix B). Singly occurring species were removed from the clusters and included in the synoptic table section labeled noncharacteristic singly occurring taxa (Appendix B).

Nonmetric multidimensional scaling ordination and permutational multivariate analysis of variance analyses

To visualize similarities and differences between and among plant communities, relative to environmental variables across the chrono-toposequence, we used the nonmetric multidimensional scaling (NMS) method (Kruskal 1964). Variation in plant community composition was calculated using Sørenson’s distance measurement ($\mathit{\beta }$ =0.25) and run on “slow and thorough” in the autopilot setting within PC-ORD 6 software (McCune and Mefford 2011). The ”autopilot” setting in PC-ORD 6 performs multiple runs, choosing the solution with the lowest final stress from a real run in each dimension and testing for significance based on a 95 percent confidence interval using a Monte Carlo test (McCune and Grace 2002). In the “slow and thorough” setting used here, 250 real and 250 randomized runs were used to determine the best dimensionality for the data, selecting one best solution for each dimensionality. Correlations between plant community composition and environmental variables were visualized by both a bi-plot diagram highlighting environmental variables that were correlated to plant community separation within the ordination space ($r^2$ $\ge$ 0.3) and the three NMS axes correlating to each environmental variable at each plot via Kendall’s tau correlation coefficient. NMS axes were labeled with environmental variables with a Kendall’s tau value greater than or equal to 0.5 ($\mathit{\tau }\ge \left|0.5\right|)$. Abiotic and biotic factors correlated with plant community separation in ordination space were tested for significant effects on plant community structure and significant interaction between factors, based on a 95 percent confidence interval, using a permutational multivariate analysis of variance analysis in R (R Core Team 2021).

Species richness and percent cover

To calculate species richness across the transect, the total number of species for each plant growth form (lichens, bryophytes, forbs, graminoids, shrubs, biological soil crust, and seedless vascular plants) was averaged from each glacial deposit on rocky and fine-grained substrates. To calculate percentage cover, field cover estimates of each plant growth form from each plot were averaged from each glacial surface across the transect.

Results

The time course of succession after deglaciation

Measurements of the largest thalli of three groups of lichen taxa, taken from ten sites across both the ELIA and NEO moraines, were used to estimate times of moraine stabilization using lichenometry. The largest R. geographicum s.l. thalli measured across the ELIA moraine ranged from 5 to 13 mm in diameter (mean = 7.8 ± 1.78 mm), which provides a maximum limiting age on ELIA moraine stabilization of 36 years since deglaciation. The largest Pseudephebe spp. and Umbilicaria spp. thalli on the ELIA moraine ranged from 29 to 36 mm in diameter (mean = 25.8 ± 3.97 mm), which provides a minimum limiting age for stabilization of the ELIA moraine of 22 years since deglaciation. Thus, the estimated stabilization of the ELIA moraine took place between 22 and 36 years ago (1995–1978 CE; Figure 4). The largest R. geographicum s.l. thalli on the NEO moraine were between 20 and 34 mm (mean = 27 ± 2.37 mm), which yields a maximum limiting age of 475 years since deglaciation. The largest Pseudephebe spp. and Umbilicaria spp. thalli on the NEO moraine ranged from 87 to 105 mm (mean = 91.8 ± 3.37 mm), which yields a minimum limiting age of 238 years since deglaciation. Thus, the estimated stabilization of the NEO moraine occurred between 238 and 475 years ago (1542–1779 CE; Figure 4). The upper boundary of age estimates for the NEO moraine (475 years ago) is within 25 years of the date of stabilization estimated by Ellis (1978) of approximately 500 years ago. Succession continues to unfold in this system for at least 25,000 years on fine-grained substrates to 125,000 years on rocky substrates when species accumulation stabilizes and percent cover becomes highly stratified.

Figure 4. Five maximum thallus radius measurements taken across both the ELIA (triangles) and NEO (squares) moraines in the Grizzly Glacier cirque from three different lichen taxa: Rhizocarpon geographicum s.l. (shown in black) and Pseudephebe spp. and Umbilicaria spp. (shown in gray). Age estimates were derived from growth rates presented in (Calkin and Ellis 1980). Red symbols represent average age estimates of the five largest thalli on both moraines.

Species composition and structure of cirque pioneer communities relative to other arctic moraine communities

We used a cluster analysis to identify distinct plant communities along the successional gradients. The seven plant community types identified by cluster analysis are represented by distinct plant taxa that are sometimes unique to only one community type, especially in early-successional communities growing on fine-grained substrates (cluster 1, Figure 5) and in ITKII foothill moraine communities (cluster 7, Figure 5). Plant community clusters within the dendrogram (Figure 5) fall into three main subgroups of community types, defined first by time since deglaciation, with early-successional communities (<500 years since deglaciation) separating clearly from mid-late-successional communities (>500 years since deglaciation), explaining about 25 percent of all variance between plant communities (Figure 5). Plant communities are then defined by substrate type (rocky/fine-grained) and elevation (alpine/foothills, Figure 5). Some plant community structural overlap does occur between substrate type (cluster 7, Figure 5) and elevation (cluster 5, Figure 5).

Figure 5. Cluster analysis of plots (n =42) from 1,700 to 800 m.a.s.l. and from 22 to 125,000 years since deglaciation in the central Brooks Range, north into the Arctic foothills. Study plot numbers are along the left axis (red = rocky plots; green = fine-grained plots). Seven main clusters of plots are defined by large red numbers. Groups with highest percentage floristic similarity are located closest to the left edge of the diagram; dissimilarity increases as subgroup nodes expand away from branch tips. Scale bar (top) shows percentage similarity for each cluster.

Based on the seven plant community types identified in the cluster analysis, we defined species composition of each community using the synoptic table (Appendix B) and identified characteristic species via indicator species analysis. Combined, the results of these analytical tools allowed us to construct a successional timeline and identify patterns and differences among plant communities occurring along the chrono-toposequence. For full descriptions of plant community type per the Braun-Blanquet method, see Appendix A.

22–500 years since deglaciation, plant community clusters 1 and 2

Primary succession begins after twenty-two to thirty-six years since deglaciation (Figure 4) upon moraine stabilization. On fine-grained sites, the forb Chamerion latifolium and dwarf shrub Salix arctica are common and uniquely confined to pioneering communities (Figure 6A, Table 2). Both of these plants utilize aerial seed dispersal, giving them an advantage in accessing freshly exposed sites, and both species are well-known colonizers of recently deglaciated forelands elsewhere in the Arctic (Mori, Uchida, and Kanda 2013). In the Grizzly Glacier cirque, two foliose lichens, Peltigera didactyla and P. venosa, are common components of fine-grained pioneer communities (Figure 6A, Table 2). Both lichens are symbiotic with cyanobacteria that fix atmospheric nitrogen, which gives them an advantage in colonizing nutrient-poor soils on freshly exposed surfaces (Thomson and Brehmer 1984). These common plant species define the dominant plant community type found on moraines of this age on fine-grained sediments.

Figure 6. Timeline of major changes to vegetation communities on both (A) fine-grained and (B) rock substrates across a 125,000 year chronosequence and 900 m altitudinal gradient.

Table 2. Plant community type summary table with characteristic species for each community type defined across the sequence.

Community type name	Community type age	Community type elevation and substrate	Characteristic species
Cluster 1: Chamerion latifolium and Peltigera didactyla	40–500 years (ELIA and NEO)	≥1,700 m.a.s.l., fine-grained	Foliose lichens Peltigera didactyla and Peltigera venosa, erect forb Chamerion latifolium, and dwarf shrub Salix arctica
Cluster 2: Lecanora polytropa and Candellariella vitellina	40–500 years (ELIA and NEO)	≥1,700 m.a.s.l, rock	Crustose lichens Candellariella vitellina and Lecanora polytropa and foliose lichen Umbilicaria cylindrica
Cluster 3: Salix rotundifolia and Racomitrium lanuginosum	500–10,000 years (NEO and ITKII-AP)	≥1,700 m.a.s.l., fine-grained	Fruticose lichen Dactylina ramulosa; rush Luzula confusa; erect forbs Draba macrocarpa, Saxifraga cernua, and Saxifraga razshivinii; and dwarf shrub Salix rotundifolia
Cluster 4: Betula nana and Masonhalea richardsonii	25,000–125,000 years (ITKI and SAG)	900–950 m.a.s.l., fine-grained	Fruticose lichens Cladonia bellidiflora and Cladonia sulphurina; foliose lichens Masonhalea richardsonii, Peltigera canina, and Sticta arctica; cushion forb Diapensia lapponica; erect forb Polygonum bistorta; grass Poa alpina; low shrubs Cassiope tetragona and Vaccinium vitis-idaea; and dwarf shrub Betula nana
Cluster 5: Porpidia flavocaerulescens and Cladonia squamosa	10,000 years (ITKII-AP)	1,400 m.a.s.l., rock	Crustose lichen Porpidia flavocaerulescens and Schaereria endocyanea, fruticose lichen Cladonia squamosa, foliose lichens Cetraria nigricans and Umbilicaria torrefecta, and mosses Schistidium apocarpum and Tortella tortuosa
Cluster 6: Arctoparmelia centrifuga and Ophioparma lapponicum	10,000–125,000 (ITKII-N, ITKI, SAG)	800–950 m.a.s.l., rock	Crustose lichens Ophioparma lapponica and Rhizocarpon cinereovirens and foliose lichen Arctoparmelia centrifuga and Asahinea scholanderi
Cluster 7: Mixed plant community cluster	25,000 years (ITKII-N)	800 m.a.s.l., fine-grained	Crustose lichen Ochrolechia upsaliensis; foliose lichens Hypogymnia austerodes and Lobaria linita; moss Tomenthypnum nitens; mat forb Dryas integrifolia; erect forbs Antennaria friesiana, Bistorta vivipara, Oxytropis jordalii, Pedicularis lanata; and sedge Kobresia myosuroides

‘

The Chamerion latifolium and Peltigera didactyla community type (Mucina et al. 2016; cluster 1, Figure 5 and Table 2) is a fine-grained, early-successional community type. These are alpine communities consisting of a low shrub layer, an herbaceous layer, and a cryptogam layer. Total vegetation cover is around 9 percent, dominated by lichens and bryophytes. This community type has one very highly characteristic species (phi ≥ 0.8, p < .05, mean phi ± SD = 0.25 ± 0.10), the erect forb Chamerion latifolium (100 percent frequency; Table 2, Appendixes A and B). The most abundant species in this community type is the moss Polytrichum juniperinum, the foliose lichen P. didactyla, and the forb C. latifolium. P. juniperinum is noncharacteristic to this community type because it is ubiquitous across all glacial deposits in this sequence (Appendix B).

At rocky sites, pioneer taxa are composed solely of lichen and bryophyte species (Figures 7B and 7D). R. geographicum s.l., Umbilicaria, and Pseudephebe spp. are the first plants to colonize recently deglaciated rock surfaces (Figure 6B). Additionally, the crustose lichens Lecanora polytropa and Candelariella vitellina are less common but consistently present in early-successional rocky communities as pioneer species. These lichens persist throughout the sequence on rocky substrates and define the dominant plant community type found on rocky substrates of young (twenty-two to thirty-six years since deglaciation) glacial deposits.

Figure 7. (A), (B) Average species richness and (C), (D) total percentage cover of major growth forms on (A), (C) fine-grained and (B), (D) rock substrates across chrono- and toposequences (ELIA =26–40 years since deglaciation, NEO =238–475 years since deglaciation, ITKII =10,000–25,000 years since deglaciation, ITKI =25,000–75,000 years since deglaciation, and SAG =125,000–780,000 years since deglaciation).

The Lecanora polytropa and Candellariella vitellina community type (Mucina et al. 2016; cluster 2, Figure 5 and Table 2) is a rocky substrate, early-successional community type. These are alpine communities consisting of a single cryptogam layer. Total vegetation cover is between 29 and 37 percent on average and is dominated by crustose lichens. This cluster has one very highly characteristic species (phi ≥ 0.8, p < .05, mean phi ± SD = 0.25 ± 0.10), the crustose lichen Lecanora polytropa (100 percent frequency; Table 2, Appendixes A and B). The most abundant species in this plant community are three foliose lichens in the genus Umbilicaria: U. cylindrica, U. proboscidea, and U. hyperborea. U. cylindrica is a characteristic species in the cluster and occurs in 100 percent of plant communities of this type, along with U. hyperborea, which is common on rocky substrates across the chrono-topsequence.

Pioneer species on the ELIA moraine are few in number and low in abundance. On average, rocky substrates host eight species of lichens and bryophytes with 27 percent cover (Figures 6B and 7D), and fine-grained sites host thirteeen species of forbs, shrubs, lichens, and bryophytes with only 9 percent cover (Figures 6A and 7C). Pioneer vascular plant species (e.g., the forb Chamerion latifolium and the shrub Salix arctica) generally do not persist into the later stages of succession; however, pioneer lichens and bryophytes commonly do. Pioneer species on rocky substrates (e.g., R. geographicum and U. hyperborea) are the most likely of all taxa to persist throughout the entire successional sequence and to occur at all elevations.

500–10,000 years since deglaciation, plant community clusters 3 and 5

After 500 years, fine-grained communities host increasingly diverse and abundant taxa of all growth forms, with graminoids, bryophytes, and lichens increasing in cover (Figures 6A and 7C, Table 3). By this time, small populations of graminoids, shrubs, and forbs have colonized rocky substrates, presumably in response to the accumulation of fine sediment in microtopographic depressions (Figures 6B, 7B, and 7D, Table 3). Both fine-grained and rocky plant communities continued to increase in their species richness and cover (Table 3), with fine-grained communities maintaining higher species richness and cover relative to rocky communities (Figure 7).

Table 3. Percentage change and number of plant species accumulated or lost from one glacial surface to the next across the chronosequence.

	ELIA to NEO (22–500 years) change in species richness		NEO to ITKII-AP (500–10,000 years) change in species richness		ITKII-AP to ITKII-N (10,000 years) change in species richness		ITKII-N to ITKI (10,000–25,000 years) change in species richness		ITKI to SAG (125,000+ years) change in species richness
	Alpine		Alpine		Alpine to foothills		Foothills		Foothills
	Species percentage change	Species added or lost	Species percentage change	Species added or lost	Species percentage change	Species added or lost	Species percentage change	Species added or lost	Species percentage change	Species added or lost
Rock	56	+10	38	+11	15	+5	15	+6	−25	−10
Fine-grained	44	+10	30	+10	13	+5	−6	−2	−8	−3

The Salix rotundifolia and Racomitrium lanuginosum community type (Mucina et al. 2016; cluster 3, Figure 5 and Table 2) contains all fine-grained plant communities from alpine glacial deposits (above 1,700 m.a.s.l.) older than forty years since deglaciation (NEO and ITKII-AP). This community consists of a low shrub layer, an herbaceous layer, and a cryptogam layer. Total vegetation cover is around 68 percent and is dominated by forbs and cryptogams. This cluster has four highly characteristic species (phi ≥ 0.5, p < .05; Table 2, Appendixes A and B). The most abundant species in this community are the lichens Dactylina arctica, Cladonia arbuscula, and Cetraria ericetorum, none of which are characteristic species, because they commonly occur in other community types across the sequence.

The Porpidia flavocaerulescens, Cladonia squamosa community type (Mucina et al. 2016; cluster 5, Figure 5 and Table 2) defines rocky substrate, mostly alpine plant communities on glacial deposits of this age. This community type consists of a low shrub layer, an herbaceous layer, and a cryptogam layer. Total vegetation cover is around 102 percent, with overlapping lichen thalli, and is dominated by lichens of all growth forms. There are seven highly characteristic species (phi ≥ 0.5, p < .05; Table 2, Appendixes A and B) in this community type. The most abundant species are crustose lichens Rhizocarpon geographicum and Rhizocarpon cinereovirens and foliose lichen Umbilicaria hyperborea, all common components of rocky substrates across the chrono-toposequence (Appendix B).

10,000–25,000 years since deglaciation, plant community clusters 6 and 7

The successional trajectory begins to change after 10,000 years of directional species richness and cover increase, and this is especially evident in the foothill communities of the ITKII fine-grained surfaces. The ITKII foothill moraine hosts plant communities that are unique in composition compared to all the other glacial surfaces in the sequence (cluster 7, Figure 5). This moraine is where mat-forming bryophytes such as Tomenthypnum nitens become prevalent and where forb and graminoid cover increases (Figures 6A and 7C). Additionally, the ectomycorrhizal mat-forming dwarf shrub Dryas integrifolia is consistently abundant on this moraine alone. The ITKII-N unique plant communities stand out due to potential functional competitive advantages (i.e., mat-forming and ectomycorrhizal) of species unique to this moraine, the peak in species richness on fine-grained surfaces that occurs here (Table 3), and relative stabilization of total vegetation cover from this point forward in the chrono-toposequence. These factors suggest that the ITKII surfaces are a transition point between earlier successional plant communities undergoing directional succession and postsuccessional communities that are becoming increasingly stabilized. Based on these shifts in the trajectory of community development, we consider all fine-grained communities ≥25,000 years since deglaciation (ITKI, SAG) to be “mature” communities within this system.

The unique fine-grained plant communities on the ITKII-N moraine defines the mixed transitional community type (cluster 7, Figure 5 and Table 2). It contains two-fifths of the rocky plots and all of the fine-grained plots of the ITKII-N foothill moraine, making it the only community type with a combination of both rocky and fine-grained plots. This community type consists of a low shrub layer, an herbaceous layer, and a cryptogam layer. On average, vegetation cover is >200 percent and is dominated by forbs and lichens (Figures 7A and 7C). Despite the inclusion of two substrate types in this cluster, there is one very highly characteristic species, the matt forb Dryas integrifolia (100 percent frequency, mean phi ± SD = 0.26 ± 0.10), and nine highly characteristic taxa, all of which are unique to this community type (phi ≥ 0.5, p < .05; Table 2, Appendixes A and B), making this community type the richest in unique taxa. Several of these taxa are characteristic of previously described vegetation types occurring in dry and mesic nonacidic habitats along floodplains of the Sagavanirktok River in the Arctic foothills (Walker 1985; Schickhoff et al. 2002). The most abundant species in this cluster are Rhododendron lapponicum and Dryas integrifolia; R. lapponicum is common on foothill moraine ridges in fine-grained communities.

On rocky substrates, the presence of seedless vascular plants is unique to ITKII-N and SAG communities, defining the cluster 6 community type (Figure 5, Table 2). On rocky sites in the alpine, shrubs increase in species richness and cover slightly after 10,000 years of succession, resembling foothill rocky communities in both species richness and percentage cover of shrubs throughout the remainder of the sequence (Figures 7B and 7D). As in the fine-grained communities, ITKII-N rocky communities show an increase in forb species richness and cover around this time (Figures 7B and 7D). Species richness on rocky substrates continues to increase over a longer period than on fine-grained substrates (Figure 6, Table 3). In older foothill rocky communities, bryophytes become less frequent and lichens become more diverse, including species such as Arctoparmelia centrifuga, Ophioparma lapponica, and Asahinea scholanderi, which are able to grow on both rock and soil substrates, indicating that flexibility in substrate preference might provide a competitive advantage via the ability to avoid overlapping thalli (Rogers 1990; Armstrong and Welch 2007).

The Arctoparmelia centrifuga and Ophioparma lapponicum community type (Mucina et al. 2016; cluster 6, Figure 5 and Table 2) is a rocky, late-successional stage community type spanning foothill moraines aging from 10,000 to 125,000 years since deglaciation (ITKII-N, ITKI, SAG). This community consists of a low shrub layer, an herbaceous layer, and a cryptogam layer. Total vegetation cover is around 109 percent, with overlapping lichen thalli (Figure 7D). This community is dominated by lichens of all growth forms. This cluster has four highly characteristic species (phi ≥ 0.5, p < .05; Table 2, Appendixes A and B). One moderately characteristic crustose lichen species, Lecidea lactea, is unique to this community type. The most abundant plants in this community type are A. centrifuga, R. cinereovirens, Parmelia omphalodes, and U. hyperborea, all common species on rocky substrates throughout the sequence.

A few notable altitudinal differences in plant community composition and cover are apparent between the two ITKII sites (10,000–25,000 years since deglaciation). First, seedless vascular species occur at foothill fine-grained sites but are absent in the alpine (Figures 6A and 7A). Though species richness of forbs is comparable at high and low elevations on fine-grained sites, cover increases from around 5 percent in the alpine to 50 percent in the foothills (Figure 7C). A similar pattern can be seen on rocky substrates (Figure 7D). Perhaps most striking is the increase in shrub cover between alpine and foothill ITKII surfaces; species richness increases some with decreasing altitude, and cover jumps from 1 percent in the alpine to 10 percent in the foothills, increasing or maintaining throughout the rest of chrono-toposequence (Figure 7C). This reflects the relatively low tolerance of shrubs to altitude-influenced environmental conditions (e.g., thaw season temperatures, substrate stability, soil conditions). Graminoids, lichens, and bryophytes are also more abundant on both substrates on ITKII alpine surfaces than they are in the foothills (Figures 7A and 7B). Plant communities at both high- and low-elevation ITKII surfaces share many of the same plant species, but communities at lower elevations are 15 percent more diverse on average (Figure 7, Table 3). Species that grow at both the alpine and foothill ITKII sites represent 36 percent of total ITKII species; 26 percent are unique to alpine and 37 percent are unique to foothill surfaces (total of 62 percent unique species; Table 3). Distinct elevational differences along this gradient are attributable in part to variation in seasonal temperatures and snow cover at alpine versus foothill sites. Although snow cover data were not obtainable for this study, data collected from both SNOTEL (National Water and Climate Conservation Staff) and iButton logs from 2006 to 2018 indicate that foothill mean temperatures are warmer than alpine temperatures by 1.42°C on average during the growing season and 1.03°C during the winter between 2007 and 2018. However, the summer warmth index averages 33.1 at Toolik Lake in the foothills as compared to 18.71 in the alpine between 2007 and 2018, suggesting that the number of days above 0°C at a given location has a stronger impact on plant community development than do average summer and winter temperatures.

25,000–125,000 years since deglaciation, plant community clusters 4 and 6

In later successional stages, species richness continues to decline at fine-grained sites and percentage cover continues to increase (>200 percent), although this increase is slow relative to early successional stages (Figures 6A, 7A, and 7C, Table 3). In particular, graminoids increase in cover within mature communities that are nearly equally composed of graminoids, dwarf shrubs, bryophytes, and lichens (Figure 7C). Dominating these communities are the dwarf shrubs Betula nana and Vaccinium vitis-idaea, the grass Poa alpina, and a variety of fruticose lichens (e.g., Flavocetraria cucullata, Masonhalea richardsonii, Cladonia spp.; Table 2, Appendixes A and B). Forb cover and species richness decrease markedly from transitional ITKII-N surfaces (Figures 7A and 7C).

Fine-grained surfaces on glacial deposits of this age are dominated by the Betula nana and Masonhalea richardsonii community type (Mucina et al. 2016; cluster 4, Figure 5 and Table 2). This community consists of a low shrub layer, an herbaceous layer, and a cryptogam layer. Total vegetation cover is highly stratified (235 percent). This community type is dominated by acrocarpous bryophytes, fruticose lichens, and shrubs. There are eleven highly characteristic species (phi ≥ 0.5, p < .05; Table 2, Appendixes A and B), including plants reflective of the acidic soils found in these communities. The most abundant species in this community type are the low shrubs Betula nana, Vaccinium vitis-idaea, Arctous alpine, and Salix phlebophylla and the grass Poa alpina (Appendixes A and B).

On rocky substrates, species richness begins to decrease after 125,000 years since deglaciation, whereas percentage cover stabilizes at around 100 percent (Figures 7B and 7D, Table 3). Rocky communities on the ITKI surface are still dominated by lichens but also include bryophytes, forbs, shrubs, graminoids, and seedless vascular plants, making these the most species-rich rocky substrates of the entire chrono-toposequence (Figures 6B and 7B). Rocky plant communities on glacial deposits of this age are highly similar to those found on 10,000- to 25,0000-year-old glacial deposits (ITKII-N, ITKI) and share the same community type (cluster 6, Figure 5 and Table 2).

In addition to quantifying plant community changes, we used an NMS ordination analysis to illustrate how changes in plant community structure relate to environmental variables across the chrono-toposequence. The final solution of the NMS analysis contained three dimensions after seventy-five iterations, resulting in a final stress of 10.43 and a final instability of 0.0. NMS explained 84 percent of the variation between plant communities across the sequence (Figure 8). Axis 1 explained 42 percent of variation among plant communities, correlating with differences in percentage total graminoid cover, evergreen shrub cover, latitude, erect forb cover, and elevation (Figure 8). Axis 2 explained 23 percent, correlating strongly with percentage biological soil crust cover, bare soil cover, soil pH, erect forb cover, and total lichen cover (Figure 8). Axis 3 (not shown) explained 19 percent. The majority of variation across communities is explained by axis 1 and axis 2 (65 percent), so visualization in two-dimensional space is effective here. For a table of correlation coefficients with all environmental variables by NMS axis, see Supplementary Table C.

Figure 8. NMS ordination grouped by (A) substrate and (B) time since exposure. Bi-plots (represented by arrows) within ordinations represent environmental factors with a Pearson’s correlation coefficient $r^2\ge 0.3$. Axes are labeled with variables having strong Kendall’s tau ($\mathit{\tau }>0.3)$ correlation coefficients. Groups marked “A” in (B) indicate alpine glacial deposits.

Strong correlations ($\mathit{\tau }\ge \left|0.5\right|$) between plant community structure and environmental factors were tested for statistical significance via permutational multivariant analysis of variance at a confidence interval of 95 percent. Plant community structure was most significantly affected by time since deglaciation (p = .001, $r^2=0.143$), followed by latitude (p = .001, $r^2=0.050$), substrate type (p = .012, $r^2=0.034$), soil pH (p = .015, $r^2=0.030$), and elevation (p = .033, $r^2=0.039$; Supplementary Table D). It is important to note that there are significant interactions between some factors, namely, Time since deglaciation * Latitude * Percentage bare soil (p = .003, $r^2$ = 0.035) and Time since deglaciation * Latitude (p = .005, $r^2$ = 0.038). However, here latitude is redundant with time since deglaciation, because each subsequent time step along the chronosequence was further north than the last. Similarly, bare soil decreased as succession progressed, suggesting redundancy with latitude and time since deglaciation. Variation between plant communities explained by statistically significant factors and interactions together totals 63 percent (total $r^2=0.628)$, highlighting the inherent complexity of this system, the effects that many small influencing factors can have together, the importance of interacting factors, and lack of comprehensive consideration of factors influencing plant community development along this sequence (37 percent unexplained variance).

Bi-plot diagrams within ordination space highlight the strongest linear environmental correlations with plant community separation in three-dimensional space ($r^2\ge 0.3$; Figure 8). Longer arrows within the bi-plot indicate stronger correlation with between-community separation and among-community similarity. In agreement with the cluster analysis, bi-plots show that plant communities were separated first by substrate type (rocky vs. fine-grained; Figure 8A), then by time since deglaciation, followed by elevation (Figure 8B).

Successional directionality, species turnover, and growth form turnover

Documenting plant species richness and percentage cover across the chrono-toposequence allows for visualization of successional directionality and turnover based on changes to plant community structure. A total of ninety-eight lichen, thirty-seven moss, three liverworts, two seedless vascular plant, seven grass, two rush, three sedge, thirty-five forb, and thirteen shrub species were identified across all study plots, totaling 200 species. Of these 200 species, total species richness for the entire sequence was composed of 49 percent lichens, 18.5 percent mosses, 17.5 percent forbs, 6.5 percent shurbs, 3.5 percent grasses, 1.5 percent sedges, 1.5 percent liverworts, and 1 percent seedless vascular plants and rushes. Plant species richness was greatest on fine-grained surfaces of the ITKII-N (foothill) moraine (~10,000 years since deglaciation), peaking at an average of 36 species (Figure 7A, Table 3), and on rocky substrates of the ITKI moraine (~25,000 years since deglaciation), peaking at an average of 34 species (Figure 7B, Table 3). The majority of the most common species in early successional community types (40–500 years since deglaciation) are also common components of mature plant communities in the foothills, especially on rocky substrates.

Low species turnover is indicated by the presence of all plant growth forms aside from seedless vascular plants throughout the sequence (Figure 7), increasing species richness at each time step on both substrates (Figures 7A and 7B), and the persistence of many pioneer and early-successional species into the late-successional stages (e.g., the moss Polytrichum juniperinum, fruticose lichen Stereocaulon paschale, and dwarf shrub Salix reticulata; Appendix B). Some growth forms like shrubs increase in both species richness and cover substantially at lower elevations, especially on fine-grained surfaces (Figure 7C). Other growth forms demonstrate some turnover at lower elevations; on older surfaces by a reduction in species richness, such as bryophytes on fine-grained substrates; or reduced species richness and cover as seen in forbs on both substrates (Figure 7, Table 3). The most consistent component of plant communities across all glacial deposits and substrates is lichen domination or co-domination in both species richness and percentage cover (Figure 7). In fine-grained communities on the oldest surface (SAG), there is a distinct increase in shrub, graminoid, and bryophyte cover but not necessarily species richness (Figures 7A and 7C).

Discussion

Plant succession in the central Brooks Range begins rapidly after deglaciation, with low-diversity, low-cover communities forming within 40 years and then continuing for at least 25,000 years on fine-grained surfaces (Figure 6A). Even when compared to other arctic settings, the duration of succession in the Brooks Range (25,000–125,000 years) is unusually long (Svoboda and Henry 1987; Mori et al. 2008). Because succession is so slow, any reset or upset of the successional process could stall the development of plant communities essential to habitat and ecological function for centuries or millennia. Relative to many other plant communities growing elsewhere on the planet, these Brooks Range morainal communities are exceptionally ancient and thus deserving of careful conservation. Even the most species-poor communities on the oldest glacial moraines are part of a successional sere and are not haphazardly assembled following recent ecological disturbances. Instead, these seemingly simple plant communities, in terms of species composition and structure, are part of a process begun many millennia earlier. Once these communities are lost, they are not easily replaced.

The nature of plant community development along this topo-chronosequence, coupled with the strong influence of abiotic environmental factors on community structure and formation (i.e., time since deglaciation, substrate, elevation, latitude, and soil pH), implies that these plant communities are highly susceptible to climate change. Warming can result in the expansion of shrubs into northern tundra ecosystems, effectively altering the conditions necessary for maintaining permafrost and other key ecological processes (Vowles and Björk 2019). However, transitions in growth form cover that we see between high- and low-elevation communities in this study demonstrate that under current arctic–alpine environmental conditions, shrub expansion will probably be restricted to lower elevation sites. This could allow newly deglaciated sites to act as refugia for flora and fauna affected by shrub expansion.

Lichens are the most speciose growth form across the entire topo-chronosequence, comprising 49 percent of all species. This contrasts with the situation at low-altitude, maritime arctic sites where dwarf or prostrate shrubs eventually assume dominance at fine-grained sites (Hodkinson, Coulson, and Webb 2003; G. A. Jones and Henry 2003; Mori et al. 2008). Some of the pioneering vascular species in the Grizzly Glacier cirque occur widely in a variety of glacial forelands, notably S. arctica and C. latifolium (Chapin et al. 1994; G. A. Jones and Henry 2003; Mori et al. 2008; Mori, Uchida, and Kanda 2013). However, unlike pioneering shrubs and forbs, many pioneering lichens and bryophytes continue as prevalent community members throughout succession. This is probably because of the ability of cryptogams to alter microclimate regimes and to collect sediment when soil is sparce, making conditions more favorable for the establishment of vascular species (G. A. Jones and Henry 2003; Klein and Shulski 2011; Nascimbene et al. 2017).

Less than 20 percent of pioneer species appear exclusively in the earliest stages of succession. This is especially true for the lichens, a growth form generally considered to only be important in the earliest stages of succession, along with other cryptogamic species (Chapin et al. 1994; G. A. Jones and Henry 2003; Klein and Shulski 2011; Nascimbene et al. 2017). Though lichens are recognized as key components in some arctic ecosystems (Cornelissen et al. 2001; Sancho, Allan Green, and Pintado 2007; Klein and Shulski 2011), the timescales over which they participate in these ecosystems have not been fully recognized prior to the results we present here. Gradients in abiotic factors correlate with the distributions of different plant communities across the study area (Figures 5 and 8). These correlations suggest that shrub and forb growth forms are especially limited by climate at high-elevation sites especially in early-successional stages and that some pioneering taxa are also limited by climatic factors at lower elevations, especially in later successional stages. This is consistent with the trends described in other arctic–alpine systems (Caseldine 1993; G. A. Jones and Henry 2003; Schumann, Gewolf, and Tackenberg 2016), but the persistence of many early-successional cryptogam taxa into later stages of succession across the elevation gradient and on both substrate types is novel.

The importance of site age makes it a key consideration when assessing the potential impacts of ongoing climate change and human disturbances on arctic–alpine landscapes. This is especially true where alpine climates are slowly becoming more reflective of lower elevation climates in the Arctic and arctic tundra vegetation is changing rapidly (Post et al. 2009; Raynolds et al. 2013; Ravolainen et al. 2020). The relative influence of time since deglaciation versus elevation on plant community structure along this gradient (Figures 5 and 8) suggests that elevation is less important in arctic–alpine settings like the central Brooks Range than it is in lower latitude mountain ranges (Bruun et al. 2006). Other arctic–alpine studies on primary succession in glacial forelands have also shown time since deglaciation to be the most significant factor influencing plant community formation, but these studies lacked an elevational gradient (Hodkinson, Coulson, and Webb 2003; G. A. Jones and Henry 2003). Here we show that despite significant elevation change, terrain age correlates most closely with variation between and similarity among plant communities in this sequence (Figures 5 and 8).

Plant succession on the northern flank of the Brooks Range is largely directional, meaning that species accumulation continues to increase with time, accompanied by low rates of turnover (G. A. Jones and Henry 2003). Directional succession is common in the Arctic (Svoboda and Henry 1987; G. A. Jones and Henry 2003; Mori et al. 2008; Mori, Uchida, and Kanda 2013), but the low species turnover (also known as directional nonreplacement) appears to be unique to noncoastal and high Arctic environments (Svoboda and Henry 1987; Mori et al. 2008). The gradual accumulation of plant species eventually assembles communities where percentage cover increases to >100 percent. This probably reflects the slow rates of autogenic community enrichment via an increase in microclimate regulation, moisture retention, sediment accumulation, and soil formation, all of which are known to be important in environments where spatial heterogeneity in the light supply is low (Liira et al. 2002; Ravolainen et al. 2020).

One of the hazards of using a chronosequence to provide insights into long-term processes of ecological succession is that it is largely blind to past events that may have differentially affected surfaces of differing ages (Chapin et al. 1994). The plant communities on the older glacial surfaces in the foothills inevitably experienced a number of disturbances that could not be accounted for in our study. For example, moraines >125,000 years old experienced multiple stadial and interstadial fluctuations in climate, which undoubtedly altered their vegetation and could have caused soil erosion and/or accumulation. We did our best to offset these biases by replicating plots across wide areas of each glacial deposit and by objectively choosing certain habitat types that could be easily identified and were only located on morainal ridges (i.e., rocky and fine-grained). Our findings reflect the processes of succession, in addition to the legacy of climatic and geomorphic disturbances that plant communities were exposed to after the moraines they now grow upon were first deglaciated. One way to minimize these legacy disturbance effects would be to replicate this study along other chrono-toposequences in the same region.

The findings of this study illustrate how early-successional communities play key roles in the successional legacy in arctic postglacial settings. Thus, pioneer species and early-successional communities should be focal points of conservation efforts because of the important roles they play in establishing and developing arctic vegetation communities following major disturbance. This includes monitoring cryptogam diversity and abundance in arctic communities to help manage shrub expansion and rapid arctic vegetation change associated with climate change (Cornelissen et al. 2001; Joly, Jandt, and Klein 2009). There are at least seventy alpine glacial cirques in the central Brooks Range of similar size and retreat status (Calkin and Ellis 1980) as the Grizzly Glacier. Forty of these cirque glaciers and their moraines were mapped during glacial surveys conducted in the late 1970s (Ellis 1978). There is tremendous opportunity here for extensive documentation and research on the successional process in cirque glacial forelands. The opportunity for witnessing pioneer succession unfold in these areas is available now but only for a limited window of time.

Author contributions

S. Kasanke primarily designed the study, conducted the fieldwork, analyzed data, and wrote the article as part of her master’s thesis. D. Walker assisted with study design, method development, fieldwork and logistics, data analysis, and manuscript revisions. F. S. Chapin III and D. Mann assisted with study design, method development, data interpretation, and article revisions.

Acknowledgments

The authors thank the two anonymous reviewers and editor who helped us improve this article. S.K. is thankful for the aid of a number of people who contributed to the success and quality of this research. Dr. Christopher Kasanke provided voluntary field assistance during the summers of 2017 and 2018. Both Dr. Fryday and Karen Dillman assisted in developing the study design and questions for this research. Dr. Martha Raynolds assisted us in creating our site map (Figure 1) and provided comments on previous versions of this article. Dr. Jim Ellis aided in logistical planning and lichenometric methodology and shared photographs and data from the Grizzly Glacier cirque collected during the late 1970s and early 1980s. Andrew Cortese provided voluntary aid identifying bryophytes collected from study sites. Carolyn Parker provided voluntary verification of vascular plant species collected from study sites. A special thanks is extended to Gena Wood for her amazing contribution as a volunteer field assistant during the 2017 field season. Gena provided many hours of data collection, quality control, and sample preparation and braved the Arctic Alaskan summer weather in a tent without a stipend.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15230430.2023.2178151.

Additional information

Funding

We are grateful for funding from The Ted McHenry Biology Field Research Fund, The Randy Howenstein Memorial Field Research Fund, and The Institute of Arctic Biology Summer Graduate Fellowship for supply, travel, and logistical funding. Helicopter support into the Grizzly Glacier cirque was provided by National Science Foundation Long Term Ecological Research program. The Michigan State University Herbarium Endowment Fund provided funding for travel and training in crustose lichen identification with Dr. Alan Fryday. This work was supported by the University of Alaska Fairbanks.