Range shifts in a foundation sedge potentially induce large Arctic ecosystem carbon losses and gains

Foundation species have disproportionately large impacts on ecosystem structure and function. As a result, future changes to their distribution may be important determinants of ecosystem carbon (C) cycling in a warmer world. We assessed the role of a foundation tussock sedge (Eriophorum vaginatum) as a climatically vulnerable C stock using field data, a machine learning ecological niche model, and an ensemble of terrestrial biosphere models (TBMs). Field data indicated that tussock density has decreased by ∼0.97 tussocks per m2 over the past ∼38 years on Alaska’s North Slope from ∼1981 to 2019. This declining trend is concerning because tussocks are a large Arctic C stock, which enhances soil organic layer C stocks by 6.9% on average and represents 745 Tg C across our study area. By 2100, we project that changes in tussock density may decrease the tussock C stock by 41% in regions where tussocks are currently abundant (e.g. −0.8 tussocks per m2 and −85 Tg C on the North Slope) and may increase the tussock C stock by 46% in regions where tussocks are currently scarce (e.g. +0.9 tussocks per m2 and +81 Tg C on Victoria Island). These climate-induced changes to the tussock C stock were comparable to, but sometimes opposite in sign, to vegetation C stock changes predicted by an ensemble of TBMs. Our results illustrate the important role of tussocks as a foundation species in determining future Arctic C stocks and highlight the need for better representation of this species in TBMs.


Introduction
The impact of climate change on ecosystem carbon (C) stocks will depend on the response of individual species and their relative roles in ecosystem function [1]. In this context, foundation species and their responses to climate change play an important role due to their disproportionately large impacts on ecosystem structure and function [2]. Yet, foundation species' impact on future C stocks remains uncertain due to their coarse representation in terrestrial biosphere models (TBMs) [3][4][5]. This is especially true in rapidly changing Arctic ecosystems, which have important feedbacks to global C cycling and climate, yet are often represented as a single or limited number of plant functional types (PFTs) in TBMs [4][5][6]. Here we investigated the influence of an important foundation species (Eriophorum vaginatum L., Cyperaceae; tussock cottongrass) on Arctic C stocks in the past, present, and future.
Tussock cottongrass is an important foundation species with a unique growth form and a pan-Arctic distribution [7]. Tussock cottongrass accounts for up to one-third of primary productivity in moist acidic tundra (∼7% of the arctic tundra biome) and is often present at lower density throughout the biome [8]. Tussock cottongrass also allocates a large proportion of its biomass belowground with belowground to aboveground biomass ratios that are 3-7 times higher than other tundra species [9] (figure S1 available online at stacks.iop.org/ERL/17/045024/mmedia). Because of its large allocation to belowground biomass coupled with limited decomposition in the cold environment tussock cottongrass forms root necromass mounds that allow it to escape the saturated and often anaerobic soils of tundra ecosystems (figure 1(a)) [9][10][11]. Increased arctic temperatures will enhance the decomposition of tussock necromass and lead to soil drying that may jeopardize the persistence and size of tussock C stocks. However, our ability to assess this claim is limited by a poor understanding of the sensitivity of tussocks to climate change as well as a poor representation of this species and its unique C storage characteristics in TBMs.
Tussock cottongrass's large belowground contribution to ecosystem C stocks is not explicitly represented in current TBMs. TBMs currently characterize tundra vegetation as a single or a limited number of PFTs (⩽2 PFTs, table S1) [4,5,12,13]. These PFTs utilize average traits that are likely not representative of E. vaginatum's unique growth form and highly productive root system [4,9,11,14]. Eriophorum vaginatum tussocks and other foundation species impact several complex, non-linear ecosystem processes (i.e. soil organic C dynamics, temporal vegetation dynamics, and plant mortality), that are known to be large sources of uncertainty in current TBMs [15][16][17]. TBMs risk non-linear predictive biases by parameterizing PFTs using average traits when species significantly differ from this average [5,12,13,18]. As a result, tussock cottongrass and other foundation species have the potential to be important for future C cycling projections if tussocks are vulnerable to rapid climate change.
Numerous lines of evidence indicate that the tussock C stock is vulnerable to future climate change. In addition to the environmental changes that disadvantage tussock formation in a warmer climate, climate change has shifted optimal environmental conditions for tussock cottongrass populations northward [19]. Field experiments indicate that both fertilization and warming cause declines in tussock density through increased competition with taller statured shrubs [11,[19][20][21]. We build upon these lines of evidence and hypothesize that recent climate change has altered tussock density and that climate-induced shifts in tussock density have the potential to significantly impact regional changes in ecosystem C stocks over the next century. Addressing these hypotheses will provide impetus to further understand foundation species and represent them in TBMs to improve future C cycle projections.

Historical data and resurveys
We resurveyed ∼38 year-old historical plots in 2018/2019 to determine whether tussock density has responded to recent climate change. Historical tussock density surveys at 20 sites on the North Slope of Alaska were collected between 1977 and 1982 (figure 1(b)) by [22,23] (table S2). The sites were resurveyed in 2018/2019 with the aid of one of the original authors (Ned Fetcher). Fetcher et al's [22] sites were permanently staked allowing for easy relocation. Shaver et al's [23] sites were relocated by digitizing historical maps [24]. We calculated the change (∆) in tussock density by subtracting our modern (i.e. 2018-2019) and historical (i.e. 1977-1982) plot level surveys (∆ = Modern-Historical). Hence, a negative ∆ indicated a decrease in tussock density and a positive ∆ indicated an increase in tussock density since the 1980s. We tested for changes in tussock density between the historical and modern surveys using a paired Wilcoxon signed-rank test because the data violated the normality assumptions for a parametric test.

Contemporaneous surveys of tussock density and soil organic C stocks
We quantified the contribution of tussocks to ecosystem C stocks with field surveys distributed across the Arctic. Tussock density was measured in 25-50 2 × 2 m quadrats placed at independent randomly selected distances along 100-200 m transects at 98 sites in the Alaskan, Canadian, and Russian Arctic (2015-2019, figure 1(c) and table S3). We randomly measured tussock diameters with tree callipers at 78% of the sites. We used diameter and tussock density measurements to calculate the tussock C stock with allometric equations. The tussock C stock included both above and belowground biomass and necromass in tussocks. The tussock allometric equation predicted the tussock C stock as a function of tussock diameter and was developed by harvesting and dissecting 57 tussocks in Russia and Alaska. The harvested tussock material was oven-dried at 60 • C for 48 h, weighed, and the C mass was quantified using a conversion factor of 0.5 gC g −1 . Tussock allometries from Russia and Alaska were statistically indistinguishable (F = 1.3; P = 0.28) and had high predictive power (R 2 = 0.86, figure S2 and table S4) [25].
The inter-tussock soil organic layer C stock was quantified at 47 sites to determine the relative contribution of tussock C to the total soil organic layer C stock. At each site, the inter-tussock soil organic layer C stock (g m −2 ) was quantified using the thickness of each soil organic layer (m), the soil organic layer's bulk density (g m −3 ), and the soil organic layer's C content (gC g −1 ). The soil organic layer was divided into as many as two layers by color and texture: Oi/Oe or fibrous organic horizon, and Oa/A or organic horizon. At each site inter-tussock soil organic layer thickness was measured with a ruler in at least five pits. Soil organic layer bulk density was measured using five soil cores that encompassed all soil organic layers. Core dimensions (m) were measured using calipers. Samples were oven-dried at 60 • C for 48 h with rocks and large roots (>2 mm) removed manually. The bulk density of each layer was determined based upon the core volume (cm −3 ) and dry weight (g). The C content (gC g dry soil −1 ) of a subset of the soil samples was determined using elemental analysis or loss on ignition and used to calculate bulk C density (gC cm −3 ; figures S3, S4 and table S5) [25]. We also quantified tussock bulk density (g m −3 ) and C content (gC g −1 ) in 58 circular cores sampled through the center of tussocks on the North Slope of Alaska. We removed other species' roots from the tussock root necromass cores based on color and morphology. Significant differences between the observed bulk C density (gC m −3 ) of tussock root necromass and the two types of inter-tussock soil horizons for the samples pooled across all sites were determined using t-tests.
We assessed whether tussocks enhance the soil organic layer C stock by comparing calculations of each site's total soil organic layer C stock using two different methods. The first method accounted for the tussock C stock whereas the second assumed that no tussocks were present. In the first method, we calculated each site's inter-tussock soil organic layer C stock (gC m −2 ) based upon our field surveys. We excluded the belowground volume occupied by tussocks (m 3 m −2 ) by assuming that tussocks are cylinders, with diameters equal to each site's average tussock diameter, that extend through the entire depth of the organic layer [10]. The total soil organic layer C stock was then calculated as the sum of the intertussock soil organic layer C stock and the tussock C stock. In the second method, we excluded the tussock C stock and instead assumed the belowground volume occupied by tussocks was replaced with intertussock soil. The enhancement of the soil organic layer C stock by tussocks is the difference between the first calculation, which includes tussocks, and the second calculation, which excludes tussocks and fills their belowground volume with soil.

Tussock and shrub C stock comparisons
We compared tussock and shrub C stocks across our sites to determine the importance of the tussock C stock relative to the shrub C stock. We did this because climate warming is anticipated to increase shrub abundance and biomass (i.e. 'shrubification') [5,26,27]. We quantified aboveground shrub C stocks at a subset of our sites that spanned our latitudinal gradient (34 sites, table S3). At each site, we measured the basal diameters of all Betula nana and Salix spp. shrubs within 20 numbers 0.25 m 2 plots, and calculated the shrub C stock using region and species-specific allometric equations and published above-to below-ground ratios [28]. The allometric equations quantified total shrub biomass from the measurements of basal diameter and had high predictive power (R 2 between 0.99 and 0.6, figure S5 and table S6). Shrub biomass uncertainty including uncertainty in the belowground/aboveground biomass ratios was quantified using Monte Carlo methods. For each site, we calculated the 'tussock to shrub index' which indicates the number of tussocks per m 2 that hold the same amount of C as shrubs per m 2 (gC shrubs m −2 /gC tussocks m −2 tussock −1 ). Hence, a value less than one indicates that on average a single tussock per m 2 holds more C than the shrubs in that same area.

Geospatial and remotely sensed datasets
We obtained 26 high-resolution gridded climate and edaphic data sets to determine the climatic and edaphic controls on tussock density biogeography [29][30][31][32]. Our study area encompassed all Arctic regions in North America and Eastern Eurasia as indicated by the CAVM bioclimatic zones [8]. Decadal and quinquennial means were used to minimize the impact of inter-annual weather variation and focus on climatic controls. The gridded data sets included decadal means of 19 bioclimatic variables over the periods 1967-1977 and 2005-2015, edaphic properties (i.e. soil pH, bulk density, soil texture [percentage sand, silt, and clay], and slope), and the quinquennial mean of the MODIS Summer Warmth Index (SWI, the sum of monthly mean land surface temperatures greater than 0 • C, averaged from 2014 to 2019) [33,34]. Climate and edaphic data sets were bilinearly interpolated to a 250 m common grid using ArcGIS Pro (see figure S6 workflow diagram). Spatial autocorrelation from downscaling likely did not impact our results because the coarsest product's resolution (∼5 km) was finer than the average nearest neighbor distance between our surveys (∼20 km).
Future climate projections were obtained from one model within CMIP5 (i.e. NCAR CCSM4.0). This particular model performs well in replicating historical climate patterns in the study area and was used in [6] model inter-comparison (see section 2.7). The projected monthly absolute changes in temperature and relative changes in precipitation between the baseline (2005-2015) and target years (2090-2100), in each scenario, were downscaled and applied to our climate data [35]. SWI was not available from CMIP5 and was stepped forward in time using a random forest that predicts SWI as a function of the 19 bioclimatic variables. This random forest was trained using 70 000 randomly selected pixels within the study area and validated using 30 000 independent pixels (R 2 = 0.97, See figures S6 and S7).

Tussock density and tussock C stock models
A machine learning ecological niche model was fit to the modern tussock survey, climate, and edaphic datasets to evaluate abiotic controls on tussock density and the tussock C stock. Ecological niche models are statistical models which are widely used to understand the complex relationships between species distribution and environmental factors [3,36,37]. Our machine learning ecological niche model used bias-corrected extraTrees regression to model tussock density as a function of the decadal means of the 19 bioclimatic variables and edaphic properties described in the previous section [38,39]. We chose to retain all of the covariates following the recommendation of Pearson et al [3], given that variable selection using permutation importance or VSUF decreased model performance. Permutation importance (i.e. the increase in mean squared error [MSE] when permutating a predictor or groups of predictors) and partial dependence analyses on the tussock density model were used to explore the relationship between tussock density and each explanatory variable [40,41]. Predicted tussock density was converted to the tussock C stock using another biascorrected extraTrees regression with MODIS SWI as an additional driver and observed tussock C stocks as the dependent variable (see figure S8 workflow diagram). The ecological niche model was validated against independent data. First, 30% of the survey data was withheld from the parameterization to assess performance. Second, we made predictions using mean gridded bioclimatic variables characterizing the decade before the earliest historical tussock density survey (i.e. 1967-1977) with static edaphic properties to test the model against observed historical tussock density (i.e. the surveys collected in 1977 and 1982). Historical tussock density was bias-corrected using a common approach wherein the difference between predicted historical tussock density and predicted present-day tussock density was added to current day observed tussock density (i.e. the resurveys collected in 2018/2019) [42]. The spatial extent of the historical surveys and re-surveys was approximately an order of magnitude lower than the transect surveys used to parameterize our model (i.e. 24 m 2 per site for [23] versus 200 m 2 per site in our surveys). This difference in spatial extent resulted in bias which is visible in a comparison between co-located historical and modern surveys (see figure S9). In other words, our 200 m surveys characterized a greater portion of the site and generally observed lower tussock density. As further validation, we directly compared raw predicted and observed changes in historical tussock density using a paired Wilcoxon signed-rank test.

Tussock density and tussock C stock projection
To quantify the importance of tussocks at the regional scale, we used the ecological niche model to predict modern and future tussock density and their C stocks. We divided the study area into regions based upon established political/ecological boundaries: the Republic of Sakha (SH), the Chukotka Autonomous Okrug (CH) the Yukon-Kuskokwim Delta (YK), the Seward Peninsula (SP), the North Slope, northern Canada (NC), and Victoria Island (VI). The ecological niche model was projected for the modern era, and into the future under RCP 4.5 and 8.5 for two time periods 2050 and 2100 using the climate drivers described above. Future ecological niche model runs were driven by climate with static edaphic conditions given the lack of datasets that project edaphic conditions. Our results are likely not sensitive to this assumption since edaphic properties are unlikely to change within the next century [3,37]. Our estimates and projections were further validated using back-of-the-envelope calculations. These calculations used predicted tussock density changes and average present-day tussock mass to ensure the estimates of our tussock C stock model fell within an expected range.

Comparison with TBMs
We contextualized our projections of future tussock C stocks by comparing them to projections made by an ensemble of TBMs. Gridded estimates of the total vegetation C pools to 2100 were obtained from five TBMs under the RCP 8.5 scenario from [6]. The model inter-comparison of [6] utilized the same scenario as our projections (RCP 8.5/ NCAR CCSM4.0). The predicted changes in vegetation C stocks between 2015 and 2100 averaged for the five models were calculated for each region from the gridded model outputs and compared to projected change in the tussock C stock derived from our ecological niche model.

Assessing ecological niche model uncertainties
Given that our modeling analysis synthesizes many streams of data to make inferences at large spatial scales, we used Monte Carlo simulations to assess and propagate uncertainty in our ecological niche model [43][44][45]. The underlying survey data was resampled with replacement and tussock C stock calculation, model fitting, and projection were permuted (n = 500). The resulting ensemble of outputs was used to calculate standard errors for our projections. The Monte Carlo simulations propagate uncertainty through our workflow including uncertainty in the surveys, upscaling of the surveys, and the two machine learning models (see figure S8 workflow diagram). This analysis allows us to assess the precision and reliability of our model projections across space and time and facilitates comparisons between our model and other data products [45]. It also allows for future work to determine if the accuracy and precision of predictions are increasing over time as new processes or methods are considered [45].

Historical tussock surveys
We found a statistically significant decline of 0.97 ± 0.31 tussocks per m 2 over the past ∼38 years across the historical sites on the North Slope of Alaska (P < 0.01) ( figure 2(a)). Changes in tussock density across the region ranged from −4.2 to +1.4 tussocks per m 2 . Tussock density significantly declined at 80% of the sites, whereas 20% exhibited slight increases ( figure 2(b)). The changes in tussock density were independent of latitude (Slope = 0.13 ± 0.56, P = 0.82, figure 2(b)), and of the sites which exhibited increases only two of four were significant at the 0.05 level. Decreases in tussock density coincided with a period of significant environmental change in Northern Alaska characterized by increases in air temperature and precipitation in the Arctic region (+0.6 • C and +1.5%-2% precipitation per decade) [46].

Tussocks contribution to near-surface C stocks
The presence of tussocks enhanced the soil organic layer C stocks by 0%-30% (mean = 6.9 ± 1.3%, figure 3(a)). The harvested tussocks were largely composed of root necromass (70% on average) and dead tiller necromass (20% on average) which locally extended the O soil horizon. The enhancement by tussocks was dependent upon tussock density, which explained 61% of the variation in the percent enhancement of the soil organic layer C stock by tussocks. Based upon this linear relationship, each additional tussock per m 2 enhanced the soil organic layer C stock by 3.9 ± 0.5% ( figure 3(a)). Tussocks enhanced soil organic carbon stocks by elevating the surface and packing almost twice as much carbon in a given volume as in the surrounding organic layer. For example, tussocks elevated the surface up to ∼10 cm and based upon our cores the bulk C density of tussocks (43.5 ± 3 kg C m −3 ) was significantly greater than that of the Oi/Oe fibrous organic horizon (21.2 ± 2.7 kg C m −3 , P < 0.01), which tussocks usually occupy at depth.
We assessed the potential for increased shrub growth to offset tussock loss by comparing the C stock of individual tussocks to shrubs in our surveys ( figure 3(b)). Based on these surveys the average 'tussock to shrub index' was 0.8 and ranged from 0.5 to 1.25 after considering uncertainty in the belowground/aboveground biomass ratio ( figure 3(b)). Therefore, at our sites, a change in tussock density between 0.5 and 1.25 tussocks per m 2 would roughly equal the quantity of C currently held in Betula nana and Salix spp. shrubs per m 2 .

Tussock biogeography
Predictions from the ecological niche model demonstrated high predictive power across space and time. The ecological niche model had a root mean squared error (RMSE) of 0.67 tussocks per m 2 (R 2 = 0.74) against independent measurements of modern tussock density and exhibited no latitudinal bias (P = 0.46, figures 4(a), S10(b) and table S7). The model had a RMSE of 1.95 tussock per m 2 (R 2 = 0.36) against historical tussock density ( figure 4(b), n = 20) and yielded an average change in tussock density that was statistically indistinguishable from in situ measurements (modeled: −0.07 ± 0.21, Historical: −0.97 ± 0.31; P = 0.07, RMSE = 1.95, n = 20).
The biogeography of tussock density was influenced by both climatic and edaphic factors. Mean temperature of the wettest quarter, pH, bulk density, annual mean temperature, and precipitation of the wettest month were the five most important predictors of tussock density in the model ( figure 4(c)). The grouped importance for bioclimatic variables (+96% MSE) was greater than for edaphic properties (+80% MSE). Temperature-related bioclimatic variables (+92% MSE) were slightly more important than precipitation-related ones (+87% MSE). The relationship between tussock density and mean temperature of the wettest quarter in the partial dependence analysis was non-monotonic with peaks around 2 • C and 10 • C (figure S11). Tussock density decreased with increasing soil pH, bulk density, and precipitation of the wettest month (figure S11). Tussock density increased up until an annual mean temperature of around −10 • C and then declined.
Extrapolating the ecological niche model to our region of interest, we found that the North Slope had the highest tussock density. This high tussock density resulted from favorable climatic and edaphic conditions in this region (2.3 ± 0.04 tussocks per m 2 , figure 5(a) and table 1). Edaphic and climate conditions on parts of the North Slope and Chukotka Autonomous Okrug represented an optimum in the biogeography of tussocks with acidic soils and a climate that was not too warm or cold to support high tussock densities. Low temperatures and less favorable edaphic conditions limited tussock density in Northern regions such as Banks/Victoria Island where tussock density was 65% lower than the North Slope and the Republic of Sakha where tussock density was 55% lower than on the North Slope. Edaphic conditions and warm temperatures limited tussock density in the Southern Regions such as the Yukon-Kuskokwim Delta where tussock density was 70% lower than the North Slope and the Seward Peninsula where tussock density was 57% lower than the North Slope.

Tussock shifts under future climate change
The ecological niche model projected substantial regional shifts in tussock density for all future climate scenarios (table 1). Changes in tussock density for the climate scenarios indicate that even moderate climate change will impact the biogeography of tussock density (table 1). We focused on the  The relationship between the percent enhancement of the soil organic layer carbon (C) stock by tussocks (i.e. the percent increase in C in total soil organic layer C stock due to the presence of tussocks) and tussock density. The dashed line represents a linear regression (n = 47, R 2 = 0.61, P < 0.001, y = 3.9 ± 0.46x-0.2 ± 1.2). (b) Histogram of the 'tussock to shrub index' (i.e. the average number of tussocks per m 2 required to offset shrub C per m 2 ) considering uncertainty in the belowground/aboveground biomass ratio of shrubs. A value above one indicates that on average a single tussock per m 2 holds less C than the Betula nana and Salix spp. shrubs in that same area, whereas a value below one indicates that a single tussock per m 2 holds more C than Betula nana and Salix spp. shrubs in that same area. most extreme scenario (2100 under RCP 8.5) to simplify reporting of results and because it is a common climate scenario for future Arctic C cycling. In 2100 under RCP 8.5, optimal climate conditions for high tussock density shifted northward with the greatest potential losses occurring on the North Slope (−0.8 ± 0.05 tussocks per m 2 ), northern Canada (−0.2 ± 0.04 tussocks per m 2 ) and the Chukotka Autonomous Okrug (−0.5 ± 0.05 tussocks per m 2 ; figures 5(b), (c) and table 1). The model predicted gains in tussock density in northern areas where tussock density is currently low, including Victoria Island (+0.9 ± 0.05 tussocks per m 2 ), and the Republic of Sakha (+0.7 ± 0.04 tussocks per m 2 ; figures 5(b), (c) and table 1).
Regional shifts in tussock density translated into potentially large changes in the tussock C stock in 2100 under the RCP 8.5 climate scenario (table 2). The ecological niche model explained 52% of the variation in the tussock C stock when tested against independent data and the model residuals were independent of latitude (figure S10(a), (c) and table S7). Assuming that tussock C decomposes completely, the model predicts tussock C losses of 85 ± 10 Tg C on the North Slope, 66 ± 13 Tg C in the Chukotka Autonomous Okrug, and 20 ± 5 Tg C in northern Canada. In the remaining four regions where tussock density is projected to increase by 2100, the model predicts tussock C stock increases of 2 ± 1 Tg C on the Seward Peninsula, 33 ± 10 Tg C in the Sakha Republic, 20 ± 3 Tg C in the Yukon-Kuskokwim Delta and 81 ± 7 Tg C on Victoria Island.

Comparison with TBMs
We contextualized the tussock C stock model projections with comparisons to vegetation C stocks projected by an ensemble of models for the RCP 8.5 scenario in 2100. The magnitude of the changes in the tussock C stock was comparable to changes in vegetation C stocks (table 2 and figure 6). The changes in   vegetation C stocks were positive and ranged from 17 to 148 Tg C. Changes in the tussock C stock were both positive and negative and ranged from −85 to −81 Tg C. The changes in the tussock C stock were beyond the uncertainty bounds of the changes in vegetation C stocks in every region except the Sakha Republic. Tussock C stock gains were larger than vegetation C gains on Victoria Island. On the other hand, losses of tussock C were of a different sign than increases of vegetation C in the Chukotka Autonomous Okrug, the North Slope, and northern Canada. The opposing signs of the changes in tussock and vegetation C stocks suggest that shifts in tussock density may offset expected vegetation C gains in some arctic regions.

Discussion
Here we demonstrate that tussock cottongrass is an important and climatically vulnerable foundation species in arctic tundra. Compared to other arctic species, such as shrubs, tussocks stored a comparatively large amount of C and enhanced soil C storage by elevating the surface and increasing bulk C density (figures 1, S1 and S3). Repeated historical measurements demonstrated that tussock density significantly declined in Northern Alaska over the past 38 years ( figure 2). Ecological niche modeling demonstrated that the biogeography of tussock density appears to be at a tipping point with moderate climate change inducing significant changes in tussock density in the past and into the future ( figure 2 and table 2).
Historical declines in tussock density coincided with a 0.6 • C decade −1 air temperature increase across the Arctic region [46]. Another 0.9 • C decade −1 air temperature increase is anticipated by 2100 under RCP 8.5 [47]. Recently, much attention has been paid to the shrubification of the Arctic and its C cycling implications [5,26,27,48]. Here we provided the first evidence of climate-induced shifts in another important arctic species with significant impacts on the future arctic C cycle. The high climate sensitivity of tussock cotton grass demonstrated here aligns well with observations of increasingly poor performance of this species under fertilization and warming experiments [19][20][21]. Moreover, work using other machine learning approaches also predicted a decline in tussockdominated communities over the next 100 years [3]. The ecological niche model had high predictive power for the modern biogeography of tussock density against independent testing data and provided a statistically indistinguishable historical average change in tussock density compared to observations ( figure 4(b)). This provided confidence in the model's ability to determine the environmental controls on tussock density across the region. The ecological niche model demonstrated complex and nonlinear relationships between temperature and tussock density with narrow optimal temperature ranges that resulted in high climate sensitivities (figures 4 and S11). These non-linear temperature responses align with recent work demonstrating northward shifts in the climate optimum for tussock cottongrass [19,49].
Temperature can impact tussock density directly by altering the fitness of tussock cottongrass and indirectly through temperature-related changes to the soil thermal environment, thaw depth, and soil moisture [10,11,50]. Indirect impacts appear as important as the direct impacts, given that surface subsidence was observed at the historical sites with the largest decreases in tussock density. The importance of these indirect effects in regulating tussock density may explain the decreased performance of the ecological niche model in the historical reconstruction. These indirect effects were not explicitly represented in the model, and therefore, the model predicted a more conservative historical tussock density change (Model: −0.07, Observed: 0.97). Although the model underpredicted the observed change in tussock density across the North Slope, the estimate was statistically indistinguishable from the observations providing confidence in our ability to scale the model across time. Regardless, further work that incorporates direct and indirect temperature impacts is needed to constrain future tussock density changes and their impact on C stocks.
We also demonstrated that shifts in tussock density are important for future C cycle assessments. We acknowledge the potential uncertainties in extrapolating the ecological niche model across space and time, but point to several independent lines of evidence indicating our conclusions are robust. First, our projected declines in tussock density on the North Slope align with observed decreases in tussock cottongrass performance in a warmer climate [19][20][21]. Second, tussock density changes will undoubtedly impact regional C stocks; because tussocks enhance soil organic C stocks and store more C than other arctic species. For example, if we scale tussock density changes observed over the past ∼38 years to the North Slope (assuming average present-day observed mass) the loss of tussock C would represent 22%-30% of the simulated annual net C sink of the tundra biome or 67% of the C held aboveground in North Slope shrubs [43,51]. The comparison between the C stock projections of the ecological niche and the TBMs in 2100 under RCP 8.5 also provides context and highlights how tussocks may offset or enhance C stock changes in a future warmer world ( figure 6). There are remaining uncertainties, such as the liability of tussock C which may turn over at decadal to centennial timescales (see supplementary methods S1, figure  S12). Further work to represent this species in TBMs will be needed to constrain and quantify these uncertainties.
Explicitly representing tussocks in TBMs will be challenging, but as highlighted by this work, necessary to understand future C stocks. We argue that tussock's disproportionate impact on ecosystem C stocks and high climate sensitivity warrant the development of a 'tussock' PFT in TBMs. This PFT must represent tussock density shifts associated with recruitment in novel and climatically favorable northern environments as well as competition with taller statured shrubs in southern regions. Field observations of tussock cottongrass indicate poor recruitment into new environments and poor performance in shade [3,26,52,53]. Since our ecological niche model does not explicitly represent these processes, it likely underestimates tussock C losses in the south and overestimates tussock C gains in the north [3]. This PFT should also represent the liability of tussock C and tussocks' impact on the soil micro-environment. Tussocks create warmer and deeper soils with higher nutrient turnover, which could impact tussock decomposition [10,54]. Microbial, biophysical, and biochemical changes resulting from shifts in tundra species composition could also impact tussock decomposition [55,56]. Our work is an important first step and the impetus to understand tussocks as climatically vulnerable foundation species that will disproportionately impact future ecosystem C stocks.

Conclusion
We demonstrate that range shifts in a foundation tussock sedge potentially induce large regional ecosystem C losses and gains in the Arctic. Tussocks enhanced soil organic layer C stocks by elevating the soil surface and increasing soil bulk C density. Recent climate change has decreased tussock density on Alaska's NS over the past ∼40 years, in accordance with previous experimental and modeling work. Comparisons between projections from a statistical niche model of tussock C stocks and an ensemble of TBMs further highlight the importance of range shifts in this foundation species in determining future C stocks. Through this work, we argue that constraining vegetation's impact on future Arctic C cycling will require the implementation of a tussock PFT in TBMs. This PFT should represent the development and maintenance of tussock necromass, competition between tussocks and shrubs, and tussock's impact on the soil micro-environment. We encourage further investigation of the role of losses and gains in tussocks and other foundation species in determining future ecosystem function.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.
Melendez, E Niklinska, H Thomas, S Unger, C Vizza, M Williams, and N Zimov, for their assistance. This work was supported by the National Science Foundation (DEB 1556772 and 2103539 to A V R, DGE 1841556 to S RC, PLR 1418010 to N F), the University of Notre Dame, Fulbright (open study/research grant to S R C), National Geographic (Young explorer grant to S R C), the Next-Generation Ecosystem Experiments supported by the Office of Biological and Environmental Research in the Department of Energy Office of Science, the Natural Sciences and Engineering Research Council (SERFC to P M L) and the Natural Environment Research Council (NE/M016323/1 to I M). We also thank the World Climate Research Programme, the Permafrost Carbon Network, Toolik Field Station, the Arctic LTER (NSF/PLR 1637459), the North East Science Station, BP Exploration Alaska, Herschel Island-Qikiqtaruk Territorial Park, and the Inuvialuit people. The data used in this publication will be made available through the NSF Arctic Data Center.