Divergence of Arctic shrub growth associated with sea ice decline

Shrub-Ring Chronologies.

We acquired both published and unpublished deciduous shrub-ring chronologies that were distributed throughout the Arctic region [SI Appendix, Table S1, and the Dryad Digital Repository: https://doi.org/10.5061/dryad.kh1893248 (57)] and covered, if possible, the entire 40-y-long period of passive microwave satellite-based estimates of Arctic SIE (1979 to the present). In order to perform a comparable study at the biome level, our synthesis focused on two shrub genera of commonly studied (21⇓⇓⇓–25, 33, 45, 48, 54) and widespread deciduous shrubs: Betula and Salix. We analyzed both shrub-ring chronologies and shrub-ring series from individual shrubs. Building chronologies at the level of site and species afforded the opportunity to examine the quality of the data, while also providing a summary of shrub growth at the level of site and species.

All shrub-ring series (i.e., raw measurements of annual ring widths after cross-dating for each individual sample) submitted by the contributors were subjected to cross-dating quality checks. Standard dendrochronological statistical tests, including 1) the expressed population signal (58) (EPS), 2) mean sensitivity (MS), 3) first-order autocorrelation [AR(1)] and mean pairwise correlation between all cross-sections (rbar.tot) were performed using the dplR package (59, 60) in R version 3.5.1 (61) (SI Appendix, Table S1). Shrub-ring data were included in our analyses if the corresponding chronologies 1) covered the common period (1979 to 2008) and 2) had an EPS [a theoretical indicator of how well the chronology represents the population mean (58)] greater than 0.75 (SI Appendix, Table S1). Among the 32 submitted chronologies, 9 were excluded because they failed to pass our threshold for minimum chronology coherence (i.e., EPS < 0.75 for detrended series) and/or because of biased sampling design. Additionally, sites were classified as either moist or dry (SI Appendix, Table S1) according to the protocol proposed by Myers-Smith et al. (22). Moist and wet soil classes were combined into the moist soil type. The dry soil moisture type was assigned when, during the warmest month of the year, the top 2 cm of the soil were dry to the touch (more than one inspection was preferred). The moist soil moisture type was assigned when surface soils were moist throughout the growing season or standing water was present during the warmest month.

Our final dataset consisted of 23 chronologies (9 Betula spp. chronologies and 14 Salix spp. chronologies), 641 shrubs (306 Betula shrubs and 335 Salix shrubs), and 753 cross-sections (SI Appendix, Table S1 and Fig. S2). In total, 20,336 growth rings (23 chronologies, 641 shrubs) were analyzed initially (Pearson’s correlations and linear regressions) for the period 1979 to 2008. On average, 33 cross-sections (min = 15, max = 77) were included in each chronology. Although some shrubs were subjected to serial sectioning (62), only root collar–stem base cross-sections (after cross-dating) were used in chronology construction. This approach was applied because previous studies have shown that ring widths at these basal locations are most sensitive to climate (63) and because many of the studies in the synthesis only collected discs at this location. Differences in mean standardized shrub growth for chronologies and for selected climatic variables (Fig. 2 B, C, E, and F and SI Appendix, Tables S3 and S4) were compared between an early (1979 to 1993) and recent (1994 to 2008) period, i.e., by splitting the entire study period (1979 to 2008) into two equal halves of 15 y. Differences in means for these two time periods were compared using Welch’s two-sample t test.

Detrending.

All analyses were performed using standardized chronologies. Shrub-ring data were detrended using the diameter regional curve standardization (dRCS) (64, 65), an approach designed to remove the effects of stem diameter on ring width while preserving long-term trends in shrub radial growth (see below). In the first step of dRCS detrending, we analyzed the effect of stem radius/radial diameter on ring widths (SI Appendix, Figs. S8 and S9). For the majority of shrubs we observed a positive effect of stem diameter on ring width (i.e., an increase in ring widths with an increasing radius). Ring width indices (RWIs) were calculated as ratios of observed-to-expected growth and site-level indices were processed to produce signal-free chronologies (66). Site-level and responder type (Fig. 1A and SI Appendix, Fig. S1) chronologies were computed as the arithmetic mean of RWI. Mean correlation between chronologies (1979 to 2008) was 0.108 and 0.4 for increasers and decreasers, respectively.

Shrub-ring standardization was performed in CRUST (65) and initially several detrending methods were tested (i.e., modified negative exponential, splines of various rigidity, dRCS, and age regional curve standardization [aRCS]). We aimed to standardize all of the shrub-ring data using the same method to avoid potential detrending bias across species and sites. The goal of standardization was to eliminate temporal variation in ring widths that can be attributed to the age or size of the shrub stem at the time that the ring was formed, while preserving low-frequency trends in shrub radial growth that could be influenced by the sea ice decline and associated changes in climate. Standardization also normalizes ring width magnitude across species and sites, allowing for aggregation of shrub-ring chronologies. After testing a variety of common growth ring data standardization methods, we settled on dRCS as the most suitable method for both theoretical and statistical reasons. Most detrending methods are applied at the level of the individual series with the aim of removing low-frequency variation (decadal-to-century scale trends), while retaining high-frequency variation (interannual variability). The low-frequency variation that is removed in this process may be derived from both intrinsic (e.g., age, size, canopy position) and extrinsic sources (e.g., climate trends). In recent years, the RCS method has gained prominence for its potential to remove the low-frequency trends derived from intrinsic sources, while retaining low-frequency trends driven by extrinsic forces, such as trends in climate (67). Briefly, RCS detrending involves 1) aggregating individual tree- or shrub-ring series that vary strongly in age and size, 2) aligning the series by age or size to define the overall effect of age or size on the ring widths (typically a spline curve), 3) removing the empirically derived age or size effect from each series through calculation of ratios or residuals of observed-to-expected growth for an individual of a given age or size. While RCS detrending is appealing as a method to examine low-frequency trends in tree or shrub growth, a number of pitfalls have been identified (68). One of the more challenging of these pitfalls is known as the “modern sample bias,” which emerges when a single age-related detrending curve is fit to a heterogeneous sample of both fast and slow growing trees or shrubs. The fast-growing individuals tend to be younger, while the slow-growing individuals tend to be older. Errors in the fit of the single detrending curve to these contrasting groups can lead to spurious trends in the resulting tree- or shrub-ring chronology. In an effort to mitigate the modern sample bias, Briffa and Melvin (68) advocated for the use of multiple curve RCS, in which separate detrending curves are fitted to fast- and slow-growing individuals. An important constraint of multiple curve RCS is that it requires a larger sample size than single curve RCS (∼50 series/curve) (65).

Development of a shrub-ring chronology is generally much more labor intensive than development of a tree-ring chronology, due to application of serial sectioning (62) and staining of thin sections (69) versus sanding of tree increment cores. For that reason, shrub-ring chronologies are generally constructed using fewer samples than tree-ring chronologies and shrub-ring chronologies that can support multiple curve RCS are rare. The dRCS method (64, 65) is a less well-known variant of RCS, in which the detrending curve is defined by the empirical relationship between ring widths and the radius of the stem at the time the ring was formed. Replacing age with radius is a means to address the more rapid decline in ring width with age in fast-growing individuals and the more gradual decline in ring width with age in slow-growing individuals. For this reason, we argue that dRCS may be the preferred standardization method when sample sizes are too small to support multiple curve RCS and when the aim of the study is to examine low-frequency growth trends that are driven by extrinsic forces (25, 70). In our shrub-ring dataset, we generally found weak trends between ring widths and the size or the age of the stem, meaning that standardization generally had a small effect on the resulting shrub-ring chronology, as observed in other shrub-ring studies (25).

The robustness of the relationship between sea ice and shrub growth was verified using so called “double detrended” time series. For this test, dRCS ring width indices and raw Pan-Arctic SIE time series were detrended using autoregressive (AR) modeling (71). Doing so, we aimed to examine high frequency variation, where interannual variability in shrub growth was related to detrended (using the same method) seasonal SIE time series. Other detrending methods were tested, including modified negative exponential curve and cubic smoothing spline, but the former failed in retaining high-frequency variation in nonstationary time series and in the latter it was impossible to identify a single spline wave length that provided a good fit to both the SIE and shrub-ring data. The AR detrending procedure was performed in the dplR package (59, 60) in R where detrended series were represented by the residuals of an AR model divided by the mean of those residuals, yielding a series with white noise and a mean of 1 (60). By default in the dplR package, the order of the autoregressive model was selected automatically and independently to each series using Akaike information criterion (AIC). AR detrending was applied to both types of time series, i.e., dRCS shrub growth and raw SIE, using the common period from 1980 to 2008. In order to assure constant variance, all detrended time series were scaled (z-score) before the analyses.

Sea Ice Data.

Pan-Arctic and regional SIE data were obtained from the National Snow and Ice Data Center, (NSIDC) (72, 73). We acquired monthly Pan-Arctic SIE data (National Oceanic and Atmospheric Administration [NOAA]/NSIDC Climate Data Record of Passive Microwave Sea Ice Concentration, Version 2; downloaded on February 24, 2017) using the following dataset link: ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/. The monthly regional Sea Ice Index (Version 3) (73) was acquired from the following source: ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/seaice_analysis/. Data are from the NSIDC Sea Ice Index (https://nsidc.org/data/seaice_index), obtained from the Defense Meteorological Satellite Program (DMSP) series of passive microwave remote sensing instruments. The Sea Ice Index is based on the datasets near-real-time DMSP Special Sensor Microwave/Imager (SSMI/I)-Special Sensor Microwave/Imager/Sounder (SSMIS) Daily Polar Gridded (25 × 25 km) Sea Ice Concentrations (https://nsidc.org/data/nsidc-0081) and the National Aeronautics and Space Administration-produced Sea Ice Concentrations from Nimbus-7 SMMR and DMSP SSM/I Passive Microwave Data (https://nsidc.org/data/nsidc-0051). For more information on platforms, instruments, and methods used for SIE data acquisition and processing please refer to Meier et al. (72) and the Sea Ice Index user guide (https://nsidc.org/data/g02135).

Regional SIE data were assigned to each shrub-ring location based on the shortest distance of each study site to open sea/ocean. Specifically, in our analyses, the following sea ice regions were used: Central Arctic, Hudson Bay, Baffin Sea, Beaufort Sea, East Siberian Sea, Kara Sea, Barents Sea, and Greenland Sea. A complete list of the regional SIE data assigned to each chronology is presented in SI Appendix, Table S2. Two gaps in monthly Pan-Arctic and regional SIE data (i.e., for December 1987 and January 1988) were filled by computing the arithmetic mean of monthly indices for adjacent years for December only (using SIE for December in years 1986 and 1988) and January only (using SIE for January in years 1987 and 1989). In order to examine more synoptic effects of sea ice extent on summer climate and shrub growth, we used seasonal SIE data instead of monthly data in both LME and SEM analyses (see below). The complete list of seasonal SIE variables used in the analyses is presented in SI Appendix, Table S9.

Finally, in order to obtain a more detailed understanding of the regional-scale relationships between sea ice and shrub growth, additional LME models and SEM analyses for each responder group were run using previous year first advance day of sea ice (pFAD, i.e., timing of regional sea ice freeze-up) and both previous and current year last retreat day of sea ice (pLRD and LRD, i.e., timing of regional sea ice retreat). Both measures were expressed as day of the year, thus were z-scored before the analyses for all sea ice regions. Sea ice data sources and equations for FAD and LRD for each region followed Stroeve et al. (74). Two sea ice concentration (SIC) thresholds (15% and 50%) were used for computing mean timing (i.e., day of the year) of sea ice retreat and advance in annual resolution for each region. Mean FAD and LRD values were obtained by taking an area-weighted average of grid cells that experienced both retreat and advance of sea ice for every year from 1979 to 2017. The regional division of Pan-Arctic sea ice was applied after Stroeve et al. (74) and assignment of each shrub-ring chronology to a certain sea ice region was similar as it was for regional SIE (SI Appendix, Table S2). LRD and FAD at the 50% SIC threshold showed greater intraseasonal variation and were better correlated with the shrub-ring data than when computed at 15% SIC. Therefore, LRD and FAD calculated at the 50% SIC threshold were retained in the final analyses. Only standardized (z-score) sea ice-related variables were used in SEM and LME analyses (see below).

Climate Data.

Mean monthly air temperature and precipitation data were acquired from the Climatic Research Unit (CRU) time series (TS) version 4.01 using 0.5 × 0.5° grids (75). We recognize that the CRU TS dataset is not perfect and that its quality decreases with increasing latitude, particularly for precipitation (76). However, in most cases, long-term instrumental data are not available and we felt that using CRU TS data for all sites would be least likely to introduce bias into our analyses. Summer temperature and precipitation (SI Appendix, Table S4) were calculated using June-to-August (JJA) mean monthly CRU TS dataset. SPEI was calculated using CRU TS 4.01 temperature and precipitation data in the spei package (v1.7) (77, 78) in R. SPEI is a multiscalar drought index that is related to the balance between precipitation and potential evapotranspiration (PET) and is commonly used in tree-ring studies. In order to evaluate drought conditions across longer periods, the monthly SPEI was calculated at shorter (1 mo) and longer (12 and 24 mo) timescales and tested in SEMs. In our final models, summer SPEI represented by the arithmetic mean of monthly June-to-August SPEI calculated at a 12-mo timescale was used. PET was calculated using the Thornthwaite equation (79).

Climate–Growth Relationships.

We began our analyses at the coarsest scale, relating Pan-Arctic SIE with shrub-ring chronologies and progressively analyzed the data at increasingly finer resolutions down to the relationships between the timing of sea ice retreat at the regional scale with the detrended ring width series of individual shrubs. We aimed to first provide a summarized understanding of shrub growth at the site-by-species level using chronologies, and to then study individual shrub responses to sea ice using linear mixed effects models. Specifically, in order to examine the potential link between sea ice and shrub growth, we first computed Pearson’s correlations between the 23 detrended (dRCS) shrub chronologies and monthly Pan-Arctic and regional SIE (SI Appendix, Tables S5 and S7) over the period 1979 to 2008. Additionally, we ran similar analyses using the entire chronology timespan (i.e., from year 1979 when the first SIE measurements were performed until the most recent year of each chronology) (SI Appendix, Table S6). The analyses were performed using monthly (from previous June to current August) and seasonal SIE data first for the entire Arctic Ocean and for the regional SIE relevant to each chronology location (SI Appendix, Table S2). Correlation P values were adjusted for multiple comparisons at the chronology level using the false discovery rate correction (80). However, shrub and chronology aggregation to each responder group (see below) was not restricted to the outcome of the correction test. Additionally, linear regression was used to assess direct relationships between each chronology and both seasonal Pan-Arctic and seasonal regional SIE (SI Appendix, Table S8). Based on these results, each chronology that was significantly correlated with either Pan-Arctic or regional SIE was aggregated into a specific responder group: increasers, chronologies that were negatively correlated with at least one monthly or seasonal SIE variable; and decreasers, chronologies that were positively correlated with at least one monthly or seasonal SIE variable (Fig. 1B and SI Appendix, Table S1). Chronologies that were not significantly correlated with at least one SIE variable were tagged as neutral (SI Appendix, Table S1) and were not included in further analysis. To statistically corroborate our grouping of shrub-ring chronologies into increasers and decreasers, we conducted a k-means cluster analysis, which is an unsupervised machine learning method of clustering a dataset into groups (81). This analysis was done using the correlation coefficients for each shrub-ring chronology with monthly Pan-Arctic SIE between previous June and current August (SI Appendix, Tables S5 and S6). The optimal number of groups was tested by minimizing within group sum of squares. This analysis confirmed the identities of the increasers and decreasers, while also revealing that the dataset supports two primary groups (SI Appendix, Fig. S3).

All shrub-ring chronologies, except for one (SI Appendix, Tables S1 and S5–S8), were significantly correlated with Arctic-wide or regional SIE. Thirteen chronologies (59% of shrubs: 201 Betula and 162 Salix) were negatively correlated with at least one monthly or seasonal SIE variable (increasers) and nine shrub chronologies (41% of shrubs: 105 Betula and 146 Salix) were positively correlated with at least one monthly or seasonal SIE variable (decreasers) (SI Appendix, Table S1). In order to identify very large (Pan-Arctic)-scale growth signals, all standardized shrub-ring chronologies that were assigned either to increasers or decreasers were averaged to create a mean responder-type chronology. Mean correlation (r.bar.tot) between all shrubs’ growth series per responder group was 0.094 (n = 363 shrubs) and 0.108 (n = 251 shrubs) for increasers and decreasers, with 0.971 and 0.967 EPS accordingly for the common period 1979 to 2008. Mean shrub age for increasers was 45.1 y (SD = 18.3), while the median was 44 y. Mean shrub age for decreasers was 39.4 y (SD = 15.8), with a median of 35 y. Thus, our threshold for young versus old shrubs was set at 40 y (Fig. 2 A and D). For each responder-type chronology (Fig. 1A), correlation coefficients over time between seasonal Pan-Arctic SIE and shrub growth were assessed using a moving correlation analyses for the common period 1980 to 2008 (Fig. 3) using the treeclim package (82) in R. These analyses allowed us to examine potential changes (both in strength and sign) of correlations between sea ice and responder-type chronologies over time. The analyses were computed using a fixed window of 20 y with 1-y moving windows. Significance of the correlation coefficients was assessed with bootstrapped resampling (1,000 iterations).

Potential differences in growth trends and climate–growth relationships (i.e., summer temperature and summer Pan-Arctic SIE) between dry vs. moist sites and Salix vs. Betula chronologies were investigated using linear regression. The differences were examined with ANOVA and slopes comparison using the “lstrends” function in the lsmeans package (83) in R. The ratio of chronologies from sites with soil moisture classified as dry versus moist was 9:4 and 4:5, respectively, for increasers and decreasers (SI Appendix, Table S1). Among increasers there were six Betula and seven Salix chronologies, whereas for decreasers there were three Betula and six Salix chronologies.

Structural Equation Modeling.

In order to examine pathways by which SIE potentially affects shrub growth and to differentiate between direct and indirect effects of SIE and climate (temperature, precipitation, SPEI) on shrub growth, we used piecewise structural equation models (35, 84) (SEMs) (Fig. 6). In the first step of SEM analyses, all shrub chronologies that demonstrated a significant correlation with SIE (i.e., 22 out of 23 studied chronologies) were pooled into two groups of responders, i.e., increasers and decreasers. For each group of responders, shrub-ring chronologies at the level of site and species were used in a separate SEM analyses. Following standard dendrochronological protocol, we used chronologies rather than individual shrub-ring series in order to enhance the common signal for a given species at a given site. Thus, the final dataset consisted of 13 chronologies (i.e., 390 mean annual RWI records) for increasers and 9 chronologies (i.e., 270 mean annual RWI records) for decreasers analyzed for the period 1979 to 2008. Due to the lack of previous year sea ice data for shrub growth in 1979, all SEM analyses using sea ice data from the year preceding growth were performed for the period 1980 to 2008.

In this first step, a conceptual model was constructed in order to investigate complex (i.e., direct and indirect) relationships among SIE, climate, and shrub growth of both groups of responders. Considering regional climate and using a priori knowledge (22), we included in our SEM summer climate variables only, i.e., June-to-August temperature and precipitation. Growth relationships with dormant season climate are occasionally reported and were not included in our analyses. To investigate a possible link between declining shrub growth and moisture limitation (22, 23, 25), we included June-to-August SPEI in the model. In order to explore potential mechanisms behind divergent shrub growth responses associated with changes in sea ice and climate, the same SEM structure for each group of responders was applied. Thus, all paths, including nonsignificant ones, were retained in the final SEMs for comparative purposes.

In order to account for possible lagged effects of sea ice on climate and therefore on shrub growth, we fit separate models for seasonal SIE conditions present during the dormant season (i.e., previous autumn, winter, and spring). To further explore the connections between sea ice conditions, climate, and shrub growth, we fit models using previous year FAD (pFAD) and LRD (pLRD), as well as current year LRD. Due to relatively high pairwise correlation coefficients (i.e., >0.5 or < −0.5) among sea ice-related variables (i.e., Pan-Arctic SIE, regional SIE, pFAD, and LRD, SI Appendix, Figs. S10 and S11) only one sea ice variable was included in each SEM analysis and LME model. Three main potential causal pathways through which each sea ice-related variable may influence annual shrub growth were investigated in the SEM analyses, namely via 1) summer precipitation, 2) summer temperature, and 3) summer SPEI. By partitioning covariance among sea ice and climate variables, we gained insight into the pathways through which declining SIE affects tundra shrub growth.

In order to incorporate variation among the chronologies aggregated to a particular responder group, LME models for RWI, summer SPEI, temperature, and precipitation were implemented in the SEMs with a random intercept for a site included at each level of each piecewise SEM. Using this structure allowed us to test for potential site-specific variability in climate–growth and sea ice–climate relationships. We assumed linear Gaussian relationships between variables included in the model. Prior to SEM analyses, we tested each variable for normality and generated univariate density plots. To improve predictive accuracy of our SEMs, only standardized data were used. Shrub growth data (response variable) were standardized per shrub using dRCS (see Detrending above). Other variables were standardized per site (mean zero, unit variance) using z-scores, except for SPEI which was already represented by an index. As well, only standardized coefficient estimates are reported for piecewise SEM results (Fig. 6 and SI Appendix, Tables S10 and S11).

In order to specify which sea ice variable best explained the variation in shrub growth, separate SEM models were fit using seasonal regional SIE as well as regional FAD and LRD variables (SI Appendix, Tables S10 and S11). We were primarily interested in testing relationships between regional sea ice variables and regional climate to elucidate the indirect impact of sea ice on shrub growth. Thus, the models with the smallest direct effect of sea ice on shrub growth were the focus of the analysis for each responder group. The models for Pan-Arctic SIE were computed as well, for comparative purposes. Since those models highlighted strong direct effects of seasonal SIE on shrub growth, and actual direct effects of SIE on tundra shrub growth are unlikely, we emphasize the results of SEMs utilizing regional sea ice variables (Fig. 6 A and B). The direct path between sea ice and shrub growth in our SEMs is treated as an indication of additional indirect effects of SIE on shrub growth, i.e., unexplained variance. For instance, SIE is expected to affect numerous climate variables that may influence shrub growth, but are not available in long-term gridded climate datasets. These variables include, but are not limited to, cloud cover, fog, irradiance, wind speed, precipitation type, and the atmospheric water vapor pressure deficit. We anticipate that inclusion of these variables would reduce the significance of the direct path between SIE and shrub growth, while also increasing the explanatory power of our SEMs.

We evaluated the fit of each SEM using χ², P value, Fisher’s C test, and Akaike’s information criterion corrected for small sample size (AICc) (35) (SI Appendix, Tables S10 and S11). The model was judged to have a good fit if P > 0.05, which indicates that model is consistent with the data (84). Shipley’s test of d-separation was used to investigate potentially missing paths in each piecewise SEM (35, 85). For each response variable in the SEM the level of variance explained was given using conditional R² (pseudo R²) that takes into account both fixed and random effects (86) (Fig. 6 and SI Appendix, Tables S10 and S11). LME models were fit using the nlme package (87), whereas SEMs were constructed using the piecewiseSEM package (35).

Linear Mixed Effects Models.

In addition to SEMs, the link between each sea ice variable and the growth rings of individual shrubs was verified by analyzing individual shrub-ring series hierarchically in LME models. By including random effects, LME models allowed us to account for various levels of potential nonindependence in our dataset. To account for the nested structure of shrub sampling designs across all study sites (e.g., different number of shrubs of a particular species sampled across sites) various random effects were included and tested in the LME models. Specifically, LMEs were constructed for: 1) each responder group using all shrubs (not chronologies) and all sea ice-related variables tested in SEMs as fixed effects and a random intercept for shrub nested within a site (LME1) (SI Appendix, Table S12); 2) each site using all shrubs and each predictor included in the SEMs as a fixed effect and year as a random effect (LME2) (SI Appendix, Table S13); 3) each responder group using detrended shrub growth ring series and detrended sea ice-related time series (LME3) (SI Appendix, Table S14), to focus on high frequency variation. LMEs for all shrubs at the responder group level used 11,239 growth rings (n = 363 shrubs; 201 Betula and 162 Salix) and 7,822 growth rings (n = 251 shrubs; 105 Betula and 146 Salix) for increasers and decreasers, respectively. LME models covered the common period of 1980 to 2008. This period was chosen in order to compare all shrubs’ growth across all sites in the most comprehensive way (i.e., for the period represented by the vast majority of shrubs included in our synthesis) (SI Appendix, Table S1). Having right-skewed distributions, ring width indices for each responder group were square-root transformed for final LME models (SI Appendix, Fig. S1 C and F). Results obtained for transformed and untransformed shrub growth data were qualitatively similar. Fixed effects were represented by the same sea ice-related variables used in the SEMs (i.e., Pan-Arctic and regional SIE, pFAD, pLRD, LRD) (SI Appendix, Table S9). All predictor variables were transformed to z-scores prior to the analysis (i.e., mean centered and scaled by the SD). The final random intercept structure for a responder level LME models was composed of an individual shrub nested within a site (LME1). This random intercept structure accounts for variation in means between unevenly replicated units and partitions between-site variation in shrub growth–sea ice relationships. It allowed us to assign the correct weight to each annual observation between shrub growth and sea ice among study sites and minimized the effect of (both low and high) extreme years in regional sea ice declines on relationships described. To avoid the risk of multicollinearity among predictor variables (SI Appendix, Figs. S10 and S11), we primarily investigated the univariate relationships between sea ice and shrub growth. In addition to the null model, we tested 11 LME models for the period 1980 to 2008 per each responder group. As our main interest was to investigate sea ice–shrub growth interactions, LME models were fit with the exact same sea ice variables used in SEM analyses. The results of LME models using all shrubs per site and species for each responder group (SI Appendix, Table S12) confirmed the divergent sea ice vs. shrub growth signal obtained when using chronologies in SEMs.

Additionally, individual LME models were run for each species-by-site combination separately using all shrubs, instead of a chronology at a level of a particular site (LME2) (SI Appendix, Table S13). A random intercept for year was included in these models. These LME models were fit using all variables used in SEMs, i.e., both sea ice and summer climate variables. The best sea ice model for each site was compared with the overall best model. The results of LME models using all shrubs per site and species confirmed the divergent sea ice vs. shrub growth signal obtained when using chronologies. Finally, LME models were run for each responder group using double detrended shrub growth series and detrended sea ice time series for all sea ice-related variables (LME3) (SI Appendix, Table S14 and Fig. S6 C, D, G, and H). This was done in order to account primarily for high-frequency variation in shrub and sea ice time series and to compare the interannual variability in both time series. For these models, both RWI after dRCS detrending for all shrubs and sea ice-related time series were used after AR detrending. In the final models, a random intercept for a site was included as such a structure produced a lower AIC than a nested random effects structure (i.e., shrub nested within a site).

All LMEs were fit with the R package nlme (86). Maximum likelihood was used for model fit estimations. A first-order autocorrelation structure was included in each model. Diagnostic plots were used for assessing the normality of residuals and random effects, as well as to evaluate the homogeneity of variance and the assumption of linearity. Model selection and ranking was performed using the “dredge” function in the MuMIn package (88) in R. For each LME model comparison, we ran a null model that included a random intercept only and the same random effects structure as the full models. For model selection, AICc was used, together with Akaike weights for model comparison (89) or ΔAIC (i.e., difference between best model and the corresponding null model). Additionally, the quality of LME models was quantified using marginal and conditional pseudo R². These model statistics were computed for each model using the “r.squaredGLMM” function in the MuMIn package (88) in R, taking into consideration possible constraints of these statistics (86).