Strong shrub expansion in tundra-taiga, tree infilling in taiga and stable tundra in central Chukotka (north-eastern Siberia) between 2000 and 2017

Vegetation is responding to climate change, which is especially prominent in the Arctic. Vegetation change is manifest in different ways and varies regionally, depending on the characteristics of the investigated area. Although vegetation in some Arctic areas has been thoroughly investigated, central Chukotka (NE Siberia) with its highly diverse vegetation, mountainous landscape and deciduous needle-leaf treeline remains poorly explored, despite showing strong greening in remote-sensing products. Here we quantify recent vegetation compositional changes in central Chukotka over 15 years between 2000/2001/2002 and 2016/2017. We numerically related field-derived information on foliage projective cover (percentage cover) of different plant taxa from 52 vegetation plots to remote-sensing derived (Landsat) spectral indices (Normalised Difference Vegetation Index (NDVI), Normalised Difference Water Index (NDWI) and Normalised Difference Snow Index (NDSI)) using constrained ordination. Clustering of ordination scores resulted in four land-cover classes: (1) larch closed-canopy forest, (2) forest tundra and shrub tundra, (3) graminoid tundra and (4) prostrate herb tundra and barren areas. We produced land-cover maps for early (2000, 2001 or 2002) and recent (2016 or 2017) time-slices for four focus regions along the tundra-taiga vegetation gradient. Transition from graminoid tundra to forest tundra and shrub tundra is interpreted as shrubification and amounts to 20% area increase in the tundra-taiga zone and 40% area increase in the northern taiga. Major contributors of shrubification are alder, dwarf birch and some species of the heather family. Land-cover change from the forest tundra and shrub tundra class to the larch closed-canopy forest class is interpreted as tree infilling and is notable in the northern taiga. We find almost no land-cover changes in the present treeless tundra.


Introduction
The Arctic is undergoing strong climate change, which is especially prominent in high latitudes (IPCC 2018 andOverland et al 2014). Climate change is assumed to be the main driver of vegetation change in arctic and subarctic regions (Walker et al 2006). Accordingly, drastic changes have been predicted with simulations, including an increase of graminoid, shrub and tree abundance and northward movement of vegetation transition ecotones (van der Zhang et al 2013, Kolk et al 2016). Indeed, remote-sensing data shows vegetation greening (e.g. Jia et al 2003, 2009, Bunn and Goetz 2006, Zeng et al 2011. However, greening trends derived from spectral satellite data are complex (Myers-Smith et al 2020) and only with the help of field-based approaches can aspects of vegetation change such as composition, structure and biomass be investigated as, for example, in many arctic vegetation studies (Forbes et  Regional-scale compositional vegetation-change analysis, thus, requires a direct combination of fieldbased measurements and remote-sensing data. Studies from vast unexplored areas in eastern Eurasia are required as vegetation change varies regionally (Walker et al 2017). Most of the published field data on the vegetation of arctic and subarctic regions are restricted to northern North America. Some land-cover investigation studies are available for north-western Siberia, southern and northwestern Chukotka , as well as eastern coastal Chukotka (Lin et al 2012), whereas central Chukotka remains largely unexplored. Unlike the North American boreal treeline represented by evergreen conifers, the Siberian treeline is formed by larch-a deciduous conifer which characterises the seasonal forest's carbon cycling and better tolerates ice blasting, fire, and excessive heat or cold conditions (Givnish 2002). Furthermore, larch forests typically have a low canopy cover, where understorey vegetation plays a greater role in the system's physical and biological properties (e.g. latent heat flux; Xue et al 2011). In contrast to most tundra-taiga areas investigated in eastern Siberia, central Chukotka has mountainous terrain. The treeline here is restricted by a combination of elevational and latitudinal limits. Within only 140 km, treeless tundra turns via forest tundra into taiga, covering a large gradient of vegetation communities, which have a patchy pattern. Therefore, this unique vegetation of central Chukotka should be investigated at a regional scale, not as separate communities.
Commonly, broad-scale arctic vegetation changes are extracted from optical satellite missions' data. Analyses of Advanced Very-High-Resolution Radiometer (AVHRR) data (Guay, 2014) show that greening in central Chukotka from 1982 to 2012 is particularly strong. Therefore, further research into vegetation change in this region is needed. Increase in Normalised Difference Vegetation Index (NDVI) is commonly interpreted as an increase in plant biomass, cover and abundance (Myers-Smith et al 2020). However, we aim to describe not only total cover change, but which taxa play major roles in it. Normalised Difference Water Index (NDWI) and Normalised Difference Snow Index (NDSI) of winter/spring scenes may derive more composition-specific differences than NDVI alone, considering vegetation moisture content and differentiating dense from open forests, respectively. It is an advantage if the classification scheme of vegetation types is not solely based on remote-sensing derived indices, but directly involves field data to ensure association of landcover change with vegetation composition. Such an approach, in contrast to others, allows the inference of vegetation compositional changes independent from expert knowledge of vegetation types in the classifying area. Furthermore, it connects, statistically, satellite spectral information of different vegetation types and field-based taxonomical composition of these types. A model built this way is suitable to apply to areas with a similar vegetation composition.
We aim to provide quantitative information on landscape-scale compositional vegetation changes in central Chukotka over 15 years starting in 2000. We built an ordination model coupling field-based foliage projective cover data with Landsat spectral indices and used the model scores to generate a landcover classification. This allowed us to map vegetation types across four focus areas for two time slices and derive land-cover change from the difference maps. The results were used to answer i) what kind of changes in vegetation composition have happen within the investigated 15 years in central Chukotka and ii) what is the magnitude of these changes?

Materials and methods
In this study we used two sources of data: field foliage projective cover of selected taxa and Landsat satellite data. The foliage projective cover data were used to build a statistical model, which formed the basis of vegetation classes using Landsat Indices as an input. The general scheme of data processing is presented in figure 1.

Field data collection and processing
The study area is in central Chukotka, NE Siberia, Russia (figure 1). Large parts of the region are characterised by mountain complexes. Central Chukotka has a continental type of climate: cold winters with average January temperatures of −30 • C, summers with average July temperature of +13 • C, short growing seasons (100 d yr -1 ) and low precipitation (200 mm yr -1 ; Menne et al 2012a, Menne et al 2012b).
The 16-KP-04 focus area is in the Anadyr Mountains. It represents mostly mountainous terrain with elevation up to 1500 m a.s.l. and open tundra. The area around 16-KP-01 at Ilirney Lake, represents the tundra-taiga ecotone, which is characterised by diverse vegetation: open tundra, shrublands and forest tundra. The highest elevation is 1000 m a.s.l. 16-KP-03 is characterised mostly by open and forest tundra and mountainous areas (700-1000 m a.s.l) Figure 1. Scheme of the data and methods used in the study. Transformed and topographically corrected Landsat data and field data on projective cover were used in the ordination (redundancy analysis, RDA) model. The scores of the RDA model were used in the k-means and fuzzy c-means classification. Derived hard clusters were applied to Landsat images from 2000/2001/2002 and 2016/2017 to predict land-cover classes from remote-sensing data. From these land-cover maps we derived change maps.
like 16-KP-01. 16-KP-02 belongs to the northern taiga biome with low mountains reaching an elevation of 300-500 m a.s.l. in the east. The river lowlands are covered by wetlands and uplands of forests and open land with visible fire scars.
For every field site, a detailed description of the vegetation was made (Shevtsova et al 2019). We used 52 sites (V01-V16, V18-V24, V26-V51, V54, V57-V58) from 57 investigated field-sites. Some sites were excluded due to their location in mountainous terrain being affected by shadows in the satellite images. During fieldwork we investigated vegetation in circular plots of 30 m in diameter. Within each plot, we sampled five 2 × 2 m subplots, with one subplot placed at the centre of the circle and the other four placed 7.5 m away in each of the cardinal directions. Foliage projective cover of tall (>100 cm) shrubs and trees was estimated for the entire plot, whereas foliage projective cover for all major taxa of the ground layer was estimated in each of the five 2 × 2 m subplots and averaged for each plot. To accommodate the absence of the same species at compared sites and to standardise data, we applied the Hellinger transformation (Legendre and Legendre 2012) to the foliage projective cover data.

Landsat data, pre-processing and spectral indices processing
Landsat collection 1 Level 2 imagery, processed by the United States Geological Survey (USGS) to surface reflectance with a 30-m spatial pixel ground resolution, provides the best available remote-sensing product covering long-time series at a landscapescale spatial resolution. From the Landsat archive we extracted cloud-free acquisitions for the four focus areas from an early (years 2000/2001/2002) and late time-slice 2016/2017 (appendix A, table A1). We used only peak-summer acquisitions collected from mid-July to mid-August. We found only one cloudless or partly cloudless peak-summer acquisition from 2002 and one from 2017 for the 16-KP-04 area, one each from 2001 and 2016 for 16-KP-01 and 16-KP-03 and one from 2000 and 2016 for 16-KP-02. Apart from summer acquisitions, we selected cloudless snowcovered springtime acquisitions for the corresponding year and covering the same area as the peaksummer acquisitions.
Due to sensor problems on Landsat-7 since 2003 and the fact that Landsat-8 was launched in 2013 we needed to use both Landsat-7 and Landsat-8 data. To provide high consistency between the spectral bands of different Landsat missions, we applied the standard Landsat-7 to Landsat-8 transformation from Roy et al (2016), but we still observed a sensor-related bias. Therefore, we empirically derived new transformation coefficients. We used concurrent Landsat-7 (12.08.2013 and14.08.2013) to Landsat-8 (13.08.2013 and15.08.2013) acquisitions and conducted random sampling of 500 000 points (only land, avoiding noisy shadowed areas) to build linear regression models for each spectral band we used in the processing: Where ETM + −like is the transformed Landsat-8 OLI to Landsat-7 ETM +, L8 OLI is the original Landsat-8 OLI, GREEN is the green, RED the red, NIR the near-infrared and SWIR1 the first short-wave infrared Landsat surface reflectance [0-1].
To the Landsat-7 and transformed Landsat-8 data we applied the topographical c-correction (Riano et al 2003) to compensate for differences in solar illumination due to topography. We derived solar-geometry related parameters (solar zenith angle and solar azimuth angle) from the Landsat metadata. The slope and aspect angles of slopes were calculated from the ASTER Global Digital Elevation (GDEM, Tachikawa et al 2011) with the same spatial 30 × 30 m pixel ground resolution as Landsat. Additionally, we masked clouds and water bodies in the Landsat spectral-band set using an optimised threshold in the NIR band.
These pre-processed Landsat spectral bands were used to calculate Landsat spectral indices (appendix C, table C1) from summer acquisitions.

Redundancy analysis (RDA) and classification approaches
The length of the first detrended correspondence analysis axis is equal to 2.55 standard deviation units meaning the field data in our case is distributed on a short linear gradient rather than a unimodal distribution (Borcard et al 2011). A suitable constrained ordination method for short gradients assuming linear species responses along the environmental gradients is redundancy analysis (RDA, Legendre and Legendre 2012). Using RDA, we related foliage projective cover data to the Landsat spectral indices by extraction of pixel values using the field sites' coordinates. The best parameters for the RDA model were chosen according to stepwise selection and analysis of variance.
We built up the classification based on RDA scores using two clustering approaches: k-means (hard classification, appendix D) and fuzzy c-means (soft classification; FCM; Bezdek 1981; appendix D). In both cases of clustering, compositionally similar data were arranged into groups according to similar behaviours of the explanatory Landsat data. In FCM classification we set a threshold of 0.6 and values exceeding this remained unclassified and were marked as uncertain (discussed in appendix D). The optimal number of statistically significant clusters for k-means classification was chosen according to the gap statistics method (Tibshirani et al 2001).
Based on hard classification we predicted cluster numbers for Landsat data from four focus areas for both early and late time slices. We obtained difference maps for each area reflecting changes in vegetation composition over the 15-year period. We interpreted the uncertainties that could appear in the land-cover classification in the study region using FCM results from 2016/2017. We validated the accuracy of the derived land-cover maps of 2016/2017 using field data from 2018 (appendix E).

General characteristics of the vegetation field data
We observed dense forests, formed by larch (Larix cajanderi Mayr) often with characteristic fire scars in the 16-KP-02 area (Shevtsova et al 2019). Sparse forest tundra stands are common for the 16-KP-01 and 16-KP-03 areas. The sparse larch stands are commonly accompanied by low and dwarf shrubs: Salix spp., Betula nana L., Ledum palustre L., Vaccinium spp., Empetrum nigrum L. Elevational transitions between forest tundra and open tundra are often covered by patches of dwarf pine (Pinus pumila (Pall.) Regel). Open tundra on gentle slopes at lower elevations (540-700 m a.s.l) is graminoid rich and accompanied by dwarf Salix and Vaccinium. Intermediate elevations in mountain areas of 16-KP-01 (700-900 m a.s.l.) and 16-KP-03 (700-800 m a.s.l.) as well as lower

Relating field data to Landsat spectral indices in the RDA model
Landsat spectral indices NDVI, NDWI and NDSI together explained 33% of the variance in the field vegetation data. All three of them explain a significant unique portion of the variance (table 1, figure  3). Adding further indices as explanatory variables to the RDA model did not improve the explained variance. Most of the variance (21%) of the field compositional data is explained by NDVI; NDWI and NDSI, explain around 7% and 5%, respectively. The Landsat spectral indices are distributed in the two RDA-axes space, where NDVI is negatively associated with the 1st RDA axis, NDWI positively with the 1st RDA axis and NDSI positively with the 2nd RDA axis. As axes 1 and 2 carry most information (29%) only these two RDA axes were retained for further analysis.
The larch closed-canopy forest is characterised by the highest NDVI values (0.75) and low NDSI (<0.6, figure 5). Despite the high biomass content, larch closed-canopy forests have generally poor biodiversity. The vegetation composition consists mostly of abundant L. cajanderi (20-70%) with moss or grass undergrowth (figure 6).
In the RDA model, larch closed-canopy forest is characterised by a combination of low axis-1 and axis-2 scores indicating high NDVI and low NDSI. NDSI increases along axis 2, indicating highly reflective surfaces covered by snow. High axis-2 scores in combination with low or moderate axis-1 scores indicate moderate to high NDVI, reflecting graminoid tundra or forest tundra and shrub tundra, respectively. The prostrate herb tundra and barren areas class is characterised by very high axis-1 scores originating from high NDWI and low NDVI caused by low biomass and dry barren land surface.
The validation of land-cover maps showed 66% agreement between class numbers predicted for validation sites in tundra and tundra-taiga from projective cover and from Landsat Indices (appendix E). The highest disagreement (19%) was obtained for prostrate herb tundra and barren areas that were misclassified as graminoid tundra when predicted from Landsat data.

Land-cover change between 2000 and 2017
Overall, the land cover in the treeless tundra (16-KP-04) has not significantly changed (<5%) from 2000 to 2017. In contrast, the forest tundra and shrub tundra class shows a major increase in area of up to 20%-25% in the tundra-taiga (16-KP-01, 16-KP-03). Larch closed-canopy forest expanded significantly in the northern taiga (16-KP-02) accompanied by an increase of forest and shrub tundra at higher elevations.
The northernmost area (16-KP-04) is today dominated by prostrate herb tundra and barren areas (59% in 2017) and graminoid tundra (36% in 2017, appendix F, figure F1). The change map (figure 7(a) between 2002 and 2017 indicates only a slight transition of graminoid tundra to forest and shrub tundra (4%).
In the typical tundra-taiga transition (  indicates the expansion of forest and shrub tundra by 21% replacing the graminoid tundra from 2001 to 2016. It occurs mostly in the drain valleys and gentle slopes around Lake Ilirney and smaller lakes. South of Lake Ilirney we observe a change from graminoid tundra to prostrate herb tundra and barren areas due to road construction. The southernmost location (16-KP-03) represents the tundra-taiga transition as well. In contrast to 16-KP-01, it has more prostrate herb tundra and barren areas (34% in 2001 and 29% in 2016, appendix F, figure F3). The change map (figure 7(c)) shows a transition from graminoid tundra to forest tundra and shrub tundra of 19% between 2001 and 2016. Almost no change in the prostrate herb tundra and barren areas is detected. The larch closed-canopy forest slightly expanded by 4% from 2001 to 2016.
The northern taiga 16-KP-02 area is covered by forest tundra and shrub tundra (39% in 2000 and 69% in 2016, appendix F, figure F4) and larch closedcanopy forest (5% in 2000 and 13% in 2016). The change map (figure 7(d)) from 2000 to 2016 indicates a 9% transition from forest tundra and shrub tundra to larch closed-canopy forest, as well as a 40% transition from graminoid tundra to forest tundra and shrub tundra.

Dataset limitations and optimisation
Our approach yielded useful information on land-cover change despite facing several methodological challenges. We focussed on changes of major zonal vegetation types. Accordingly, azonal vegetation types were not investigated during  fieldwork including polygonal tundra in lake depressions and floodplain vegetation, which are classified as graminoid tundra or forest and shrub tundra.
The vegetation survey plots were placed within an as homogenous as possible environment to represent broad land-cover classes and allow upscaling. However, heterogeneity in tundra and tundra-taiga landscapes is often higher than the spatial resolution of most of the remote-sensing products (Myers-Smith et al 2011, Frost et al 2014. Consequently, the spectral signal of a Landsat pixel with a spatial resolution of 30 × 30 m is mixed if different vegetation types are present within the pixel. For instance, the spectral signal of a green forest canopy with low red reflectance due to photosynthetic pigment absorption and high multiple NIR scattering and thus high NDVI values can be mixed with the high red and low NIR reflectance of a lichen understorey reducing NDVI for the plot (appendix G, figure G1). Also, Pinus pumila shrubland, a characteristic vegetation type in the investigated region, could hardly be distinguished by Landsat because it is usually found alongside sparse vegetation or bare ground and covers small patches of 25-800 m 2 . For example, dense healthy green P. pumila patches occur in V23 (16-KP-03), but a large percentage of the cover in the site is of white-coloured lichen. This causes a significant decrease in NIR and an increase in red reflectance, which results in lower NDVI values for the 30 × 30 m pixel (appendix G, figure G2).
Landsat-derived NDVI and NDWI selected for building the RDA model are known for their potential to reflect vegetation-related properties. The combination of red and NIR spectral bands used for the NDVI calculation highlights plant biomass and its chlorophyll content (Tucker 1979, Boelman, 2003. NDWI utilises green and NIR reflectance and reflects the water content (Lara, 2018). NDSI has been developed to indicate snow cover (Volovcin, 1976), although it has also been used for identifying dense forests without leaves, for example, in central Alaska (Hall et al 1998). In our data, spring-NDSI values of snow-covered landscapes helped to separate closed-canopy forests that are not covered by snow due to high standing wood from other vegetation types.
Our transformation coefficients allowed a high consistency between Landsat-7 and Landsat-8 spectral bands. A comparison of Landsat-7 data with Landsat-8 data without our transformation or with application of the standard transformation (Roy et al 2016) resulted in inaccuracy due to sensor differences of the same strength as vegetation changes. We observe less variability in the spectral bands of Landsat-7 data than Landsat-8 data due to the much lower radiometric sensitivity of the Landsat-7 sensor (USGS, 2019a(USGS, , 2019b. We significantly corrected Landsat data for shadowed areas in our classification by combining two approaches. First, we used the standard way of avoiding topographic effects and variations in the sun illumination angle by employing normalised difference spectral indices instead of single spectral band values or single band ratios (Mróz and Sobieraj 2004).
Second, we applied topographical corrections (Riano et al 2003), because in the mountainous parts of Chukotka, the quality of the spectral indices was still affected by steep slopes and very strong shadows, especially for the snow-covered acquisitions.
We validated the derived land-cover maps with new field data from the expedition in 2018 (appendix E). 66% of validation sites were correctly assigned to one of the four land-cover classes. The highest misclassification, namely systematic overestimation, we found for prostrate herb tundra and barren areas over graminoid tundra. However, the derivation of relative changes compensates for the misclassification error.

Vegetation changes from 2000/2001/2002 to 2016/2017 4.2.1. Tree infilling with no evidence of substantial treeline advance
The expansion of the sole tree species Larix cajanderi is attributed to a change from forest tundra and shrub tundra to larch closed-canopy forest. The most prominent results are obtained for the northern taiga area (16-KP-02). The larch closed-canopy forest increase resulted in an absolute L. cajanderi cover increase from 2000 to 2016 of about 5% (or an area of 15.8 km 2 ) with respect to its abundance in each of the accounted land-cover classes (appendix H, tables H1, H2). The spatial pattern of larch expansion derived from the land-cover classes' differences suggests tree infilling behaviour rather than treeline advance since the larch closed-canopy forest only significantly increased in the northern taiga where it was previously present at the landscape scale. Contrarily, we found no evidence of significant larch forest advance in the tundra-taiga ecotone.
We compared the increase of larch tree cover in our study region to other Siberian subarctic regions where the treeline is formed by Larix spp. ( figure  8). Both areas of the tundra-taiga ecotone (16-KP-01 and 16-KP-03) have similar increases to northwestern Siberia as explored by . Increase of larch cover in northern taiga (16-KP-02) is twice as high as in the tundra-taiga ecotone. Despite these similarities, one needs to be careful in such comparisons since tree cover is increasing nonlinearly and one needs to consider the period of investigation (Frost & Epstein: 1965-2009this study: 2000-2017 and the different spatial resolutions of the studies. It is likely the ecological niche of L. cajanderi will shift further north with recent warming. However, larch migration to the north might be limited by seed dispersal and reproduction rates (Wieczorek et al 2017a). Our inferences of strong infilling in southern areas and a simultaneously slow treeline advance in northern areas is in accordance with simulation results of an individual-based larch model (Kruse et al 2019a). Climate-driven infilling might be overestimated in areas recently heavily impacted by fires. However, we selected a region in the northern taiga area with no fire dynamics within the investigated period.

Shrub expansion and graminoid decline
Throughout the explored region forest tundra and shrub tundra increased in area from 13% in the early time-slice to 32% in the late time-slice (appendix H, table H1). This increase in vegetation biomass is also captured by the significant greening trend between 2002 and 2008 observed in our study area, interpreted as an increase in growing season NDVI (MODIS and AVHRR, in Guay et al 2014) and increase in peak-summer MODIS NDVI from 2000 to 2018 (appendix B).
We interpret the areal increase as expansion of the forest tundra and shrub tundra into areas formerly covered by graminoid tundra. The forest and shrub tundra class has an average shrub cover of 60%, which is much higher than the shrub cover in all other vegetation classes (figure 9). The main taxa contributing to shrubification are deciduous dwarf shrubs Vaccinium vitis-idaea, Betula nana (strong contribution only in tundra-taiga), Ledum palustre and V. uliginosum (tundra-taiga, northern taiga), Empetrum nigrum (stronger contribution in northern taiga, less strong in tundra-taiga), Cassiope tetragona and Salix spp. (minor contribution in both tundra-taiga and northern taiga) and Alnus viridis ssp. fruticosa (strong contribution only in tundra-taiga, figure 10).
The shrubification goes along with a significant decrease in the graminoid tundra class from an average of 60% in 2001%-41% in 2016 throughout the explored region (appendix H, table H1). Dominant taxa of the graminoid tundra are Poaceae and Cyperaceae, which form 37% of the cover of this class, much higher than for other classes.
Myers-Smith et al (2011) describe how shrub expansion in the tundra-taiga ecotone in recent years can be related to various environmental drivers including climate (earlier snowmelt alongside higher temperatures), natural disturbances (permafrost degradation, fire) and anthropogenic impact. Many studies report tall (mostly Alnus, rarely Pinus pumila), low and dwarf (Salix, Betula) shrub expansion all over the arctic region (Sturm et al 2001, Tape et al 2006, Myers-Smith et al 2011, Frost et al 2014. Most of these studies investigate different and/or longer time periods or focus on local-scale changes (floodplains, river slopes etc). The rates of shrubification in the central Chukotka mountainous area appear to be much lower than in the other studied circum-Arctic regions that represent lowland landscapes. We suspect that most of the shrubification in the tundra-taiga transition area in Chukotka could be related to climate changes, because anthropogenic impact is low except in the 16-KP-01 region, where  (1) larch closed-canopy forest, (2) forest tundra and shrub tundra, (3) graminoid tundra and (4) prostrate herb tundra and barren areas. Shrubs (high, low and dwarf) include: A. viridis ssp. fruticosa, P. pumila, B. nana, Salix spp., V. vitis-idaea, V. uliginosum, L. palustre, Rosa sp. and Rhododendron adamsii. road construction along Lake Ilirney is indicated as a decrease in vegetation cover ( figure 7).
Trends in the annual mean air temperature suggests an increase of 0.25 • C-1 • C per decade  years trend) in our study area, as well as a decrease in snow cover duration (1-2 d per decade, 1979-2009 years trend; Liston and Hiemstra 2011). The increase in air temperature has caused warming of permafrost in recent years (Biskaborn et al 2019), which may also support a deepening of the active layer and thus can benefit plants with long root systems. According to Fyodorov-Davydov (2008), active-layer thickness significantly increased in eastern Siberia (Northern Yakutia)) from 1999 to 2006 by about 35-70 cm in typical tundra, 7 cm in the tundra-taiga ecotone and 13 cm in the northern taiga.
Despite modern climate conditions in the tundra-taiga ecotone being favourable for both shrub and forest expansion, at first, fast reproduction rates allow low and dwarf shrubs to rapidly colonise the area. These areas will be subsequently replaced by generally slower expanding forests (cf Kharuk et al;2006, Montesano et al;, Wieczorek et al 2017b, Kruse et al 2019a. Overall, our findings of shrub expansion and decline of graminoid cover agree with modelling results. For example, a simulation experiment by Carlson et al (2018) for the subarctic claims an increase in the abundance of deciduous shrubby Salix and a decline of the graminoid Carex under a warming treatment.

Stable treeless tundra
In contrast to the fast-changing northern taiga and tundra-taiga ecotone we observed a much slower change in the tundra. Even though there is only 40 km between 16-KP-01 (tundra-taiga) and 16-KP-04 (tundra), the tundra area stayed stable during the explored 15 years. However, we could detect shrub cover increase in river valleys that suggests favourable environmental conditions for shrub growth (nutrient availability, moisture, soil development, active layer depth, shelter against harsh wind). Valley shrubification is common for many arctic regions (e.g. Epstein 2013, Naito andCairns 2014), but the rates of changes are much lower in Chukotkian tundra (16-KP-04).
Similarly, stable tundra, where shrub cover is represented by Salix glauca, Betula nana and Vaccinium uliginosum, has been reported for Disko Island in west Greenland (Callaghan et al 2011). There have been few studies of Siberian tundra and many of these focus on areas of expected strong shrubification and treeline advance. Vegetation turnover might still be occurring but is hidden within our broad vegetation classes.
The circum-Arctic NDVI change supports our findings of weak greening in the Northern Chukotka tundra region, which contrasts greatly to the strong tundra greening trend on the North Slope in Alaska (Epstein et al 2018).

Conclusions
We presented a land-cover classification and landcover change over 15-years in the recent past from four areas of central Chukotka (eastern Siberia, Russia) ranging from northern taiga to treeless tundra. In our approach, foliage projective cover from vegetation plots is statistically related to Landsat spectral indices. The derived model was then used for classification of Landsat images which enabled us to detect 15-year changes in land-cover classes and interpret the magnitude and major taxonomical contributors of vegetation change.
Our analyses yielded that the class larch closed-canopy forest significantly increased only in the northern taiga where it was previously present at the landscape scale. Contrarily, we found no evidence of significant larch forest advance in the tundrataiga ecotone. This led us to the conclusion that tree infilling behaviour rather than treeline advance characterises the study area.
We observed a strong increase of the forest tundra and shrub tundra class in northern taiga and tundrataiga. Major contributors of shrubification are alder, dwarf birch and some species of the heather family. It is concomitant with a decrease in graminoid tundra which is consistent with a previous warming treatment experiment (Carlson et al 2018).
We found that treeless tundra stayed rather stable within the investigated period, which contrasts with the strong tundra greening on the North Slope in Alaska (Epstein et al 2018), but is similar to Disko Island in west Greenland (Callaghan et al 2011).
Absolute rates of vegetation change across circum-Arctic studies are only partly comparable because of differences in methodological approaches and analysed time spans, and would require a harmonised circum-Arctic approach in the future that numerically links field and remote-sensing data.

Acknowledgments
This study has been supported by the German Federal Ministry of Education and Research (BMBF), which enabled the Russian-German research programme 'Kohlenstoff im Permafrost KoPf ' (Grant No. 03F0764A), by the Initiative and Networking Fund of the Helmholtz Association and by the ERC consolidator grant Glacial Legacy of Ulrike Herzschuh (Grant No. 772852). Birgit Heim acknowledges funding by the Helmholtz Association Climate Initiative REKLIM. We thank our Russian and German colleagues from the joint German-Russian expedition 'Keperveem 2016' for support in the field. Special thanks to Lena Popova for providing technical support with Landsat data preparation.

Data availability statement
The data that support the findings of this study are openly available via the following link https://doi.pangaea.de/10.1594/PANGAEA.908570 and cited in the references as Shevtsova et al (2019).

Appendix A. Detailed description of Landsat acquisitions
Here we describe in detail the Landsat data used for retrieving the Normalised Difference Vegetation Index (NDVI), Normalised Difference Water Index (NDWI) and Normalised Difference Snow Index (NDSI). NDVI and NDWI were retrieved from peaksummer acquisitions and NDSI from snow-covered acquisitions (table A1). The pixel size of each Landsat image is 30 × 30 m. Before indices calculation the Landsat data was topographically corrected. The subsets that we used for land-cover classification were cloud free and cloud-shadow free. Additionally, we masked all water bodies. Latdsat-8 data were transformed to Landsat-7-like (see section 1.2 Landsat data, pre-processing and spectral indices processing).
Fuzzy c-means (FCM) classification results with a threshold set to 0.6 (figure D1) gave six sites (10%) out of the 52 field sites that remained unclassified and fell into the uncertainty class. These are the sites with the most mixed vegetation. The majority of the field data (90%) is strongly associated with one of the four land-cover classes. The classes are similar to k-means results (appendix F) but given as probabilities.
The mapped FCM classification is useful to understand which portion of the land surface does not belong to one of the four classes defined by hard k-means classification but instead stays unclassified. About 30% of the land surface in the study area remains unclassified by FCM: these are (1) types not represented by field data (high mountains, wetlands) and (2) transitional open-forest tundra zones.
In the mountainous parts, uncertainties mostly occur in high-elevation areas with extreme values of spectral reflectance (16-KP-04, 16-KP-03, figure  D2). These land surfaces fall outside the established classes because no field data were sampled there. In the k-means classification (appendix F) these areas are included in the prostrate herb tundra and barren areas class.
Transition areas between grass-dominant and shrub-dominant vegetation types as well as between open and forest tundra are responsible for much of the uncertainty (16-KP-01, 16-KP-03, figure D2). In the k-means results these areas are represented by both graminoid tundra and forest and shrub tundra.
In the northern taiga (16-KP-02, figure D2) more area was classified as uncertain: mostly formerly burnt areas and wetlands. Pioneer vegetation that occupied the burned areas in the first decades comprises highly abandant foliar Alnus viridis shrubs (Lantz, 2010) and dense green understorey which are difficult to separate from other vegetation groups using remote-sensing information as it just captures high-biomass green vegetation. These areas belong to larch closed-canopy forest or forest tundra and shrub tundra in the k-means classification.  Wetlands are characteristic of river floodplains and were not surveyed in the field and therefore fall mostly under uncertain. They are mostly represented as forest tundra and shrub tundra with some patches of larch closed-canopy forest in the k-means classification.

Appendix E. Validation of land-cover maps
We validated the derived land-cover maps with new field data from the expedition 'Chukotka 2018' (Kruse et al 2019b), that revisited the regions 16-KP-01 (tundra-taiga transition, 27 sites) and 16-KP-04 (tundra, 5 sites). In 2018, we did not collect field data from northern taiga, so it was not possible to validate closed-canopy forest class with this approach. We used almost the same sampling scheme during both expeditions in 2016 and 2018. Tree and tall-shrub (>100 cm) cover was estimated the same way. The difference was only in the placement of the 2 × 2 m subplots for low-shrub and ground vegetation estimation. We used a 15 m-radius at each plot but in 2018, rather than systematic placement of the subplots in the cardinal directions and the centre, we placed three subplots into each of the differentiated three major vegetation classes (e.g. lichen-dominant communities, low-shrub tundra communities with Vaccinium spp. dominance etc). We estimated foliage projective cover for the ground layer taxa and averaged, Figure D1. Results of FCM classification on 2 RDA-axes using NDVI, NDWI and NDSI as predictors. Cluster description (represented by different colours): (1) larch closed-canopy forest, (2) forest tundra and shrub tundra, (3) graminoid tundra and (4) prostrate herb tundra and barren areas. Membership is presented as the proportion of a site that belongs to a certain class.  To validate the land-cover maps, we projected the new field data from 2018, as well as corresponding Landsat spectral indices data, into the ordination space we previously created. Predicted classes were similar for the majority of 21 plots (66%, table E1). Two sites were misclassified as larch closed-canopy forest in the foliage projective cover data due to high larch cover (up to 20%), whereas they were correctly classified into the class forest tundra and shrub tundra by Landsat spectral data. Only three plots were differently classified as graminoid tundra or forest tundra and shrub tundra and vice versa. These are characterised by mixed typical vegetation, which can hardly be distinguished by the applied classifiers. While the predicted vegetation class based on the foliage projective cover data was correctly assigned to prostate and herb tundra: the Landsat spectral indices misclassified the majority (6 of 7) of open tundra communities dominated by Dryas octopetala and Cassiope tertragona into graminoid tundra. Therefore, one needs to interpret changes between these two classes with caution. However, we did not find big changes between graminoid tundra and prostrate herb tundra in tundra-taiga or tundra zones. In addition, considering systematic overestimation of graminoid tundra in both compared years compensates for the errors of misclassification. Figure F1. Colour-coded land-cover classes of 16-KP-04 (treeless tundra): a) for 2002, where dominant classes are prostrate herb tundra and barren areas (63%) and graminoid tundra (36%); b) for 2017, where dominant classes are also prostrate herb tundra and barren areas (59%) and graminoid tundra (39%). The forest tundra and shrub tundra class appears here at a very low proportion (1%). Figure F2. Colour-coded land-cover classes of 16-KP-01 (northern border of tundra-taiga ecotone): a) for 2001, where the dominant class is graminoid tundra (92%), while prostrate herb tundra and barren areas occurred in the hilly areas (5%) and forest tundra and shrub tundra is sparsely distributed in the valleys and gentle slopes (3%); b) for 2016, where the dominant classes are graminoid tundra (69%) and forest tundra and shrub tundra (27%) distributed in the valleys and gentle slopes. The prostrate herb tundra and barren areas class (3%) is still found in the hilly areas but now also includes road construction. Figure F3. Colour-coded land-cover classes of 16-KP-03 (southern border of tundra-taiga ecotone): a) for 2001, where dominant classes are graminoid tundra (57%) and prostrate herb tundra and barren areas (34%). Forest tundra and shrub tundra occurs along the waterways and on the gentle slopes of the lakes (8%); b) for 2016, where the dominant classes are graminoid tundra (42%), forest tundra and shrub tundra (28%) and prostrate herb tundra and barren areas (29%). Larch closed-canopy forest (<1%) can barely be seen as patches along the river. Figure F4. Colour-coded land-cover classes of 16-KP-02 (northern taiga): a) for 2000, where the dominant classes are graminoid tundra (55%) and forest tundra and shrub tundra (39%). To the east one can see a mountain complex where graminoid tundra turns into prostrate herb tundra and barren areas (1%). Larch closed-canopy forest (5%) is better represented in this area than in the other explored areas; b) for 2016, where the dominant classes are forest tundra and shrub tundra (69%), graminoid tundra (17%) and larch closed-canopy forest (13%). The prostrate herb tundra and barren areas class is poorly represented (<1%). Figure G1. At site V26 (16-KP-03) the lichen understorey significantly lowers the most important vegetation index (NDVI) compared with other tree-covered areas with grass or shrub understorey. While having high percentages of green healthy tree cover, the understorey is mostly covered with highly red-reflecting lichens (and green-reflecting Larix cajanderi needles) (360x180 degree panoramic image, Stefan Kruse).