Thermo-erosion gullies boost the transition from wet to mesic tundra vegetation

Continuous permafrost zones with well-developed polygonal ice-wedge networks are particularly vulnerable to climate change. Thermo-mechanical erosion can initiate the development of gullies that lead to substantial drainage of adjacent wet habitats. How vegetation responds to this particular disturbance is currently unknown but has the potential to significantly disrupt function and structure of Arctic ecosystems. Focusing on three major gullies of Bylot Island, Nunavut, we estimated the impacts of thermoerosion processes on plant community changes. We explored over 2 years the influence of environmental factors on plant species richness, abundance and biomass in 62 low-centered wet polygons, 87 low-centered disturbed polygons and 48 mesic environment sites. Gullying decreased soil moisture by 40 % and thaw-front depth by 10 cm in the center of breached polygons within less than 5 years after the inception of ice wedge degradation, entailing a gradual yet marked vegetation shift from wet to mesic plant communities within 5 to 10 years. This transition was accompanied by a five times decrease in graminoid above-ground biomass. Soil moisture and thaw-front depth changed almost immediately following gullying initiation as they were of similar magnitude between older (> 5 years) and recently (< 5 years) disturbed polygons. In contrast, there was a lag-time in vegetation response to the altered physical environment with plant species richness and biomass differing between the two types of disturbed polygons. To date (10 years after disturbance), the stable state of the mesic environment cover has not been fully reached yet. Our results illustrate that wetlands are highly vulnerable to thermo-erosion processes, which drive landscape transformation on a relative short period of time for High Arctic perennial plant communities (5 to 10 years). Such succession towards mesic plant communities can have substantial consequences on the food availability for herbivores and carbon emissions of Arctic ecosystems.


Introduction
Warming in the Arctic is occurring twice as fast as the global average (USGCRP, 2009;New et al., 2011;NOAA, 2014).Perennially frozen ground (permafrost) has consequently warmed by 2 • C over the last 20 to 30 years (Christiansen et al., 2010;Romanovsky et al., 2010), and there is now evidence of a decrease in both permafrost area extent across the Northern Hemisphere and permafrost thickness at the local scale (Beilman and Robinson, 2003;Payette et al., 2004;Smith, 2011).
Permafrost is tightly associated with biophysical components such as air temperatures, soil conditions, surface water, groundwater, snow cover and vegetation (Jorgenson et al., 2010;Sjöberg, 2015).Permafrost impedes water to drain to deeper soil layers and maintains a perched water table and saturated soils which favors the existence of wetlands (Woo, 2012;Natali et al., 2015).Permafrost degradation that would increase subsurface drainage and reduce the extent of lakes and wetlands at high latitudes (Avis et al., 2011;Jorgenson et al., 2013;Beck et al., 2015) would thus have major consequences on ecosystem structure and function (Collins et al., 2013;Jorgenson et al., 2013).It would also strongly influ-Published by Copernicus Publications on behalf of the European Geosciences Union.
Several forms of ground and massive ice can be found within permafrost (Rowland et al., 2010), especially ice wedges in regions where winter temperatures enable thermal contraction cracking (Fortier and Allard, 2005;Kokelj et al., 2014; M. T. Jorgenson et al., 2015;Sarrazin and Allard, 2015).Continuous permafrost zones with well-developed polygonal ice-wedge networks are particularly vulnerable to climate change because ice wedges are usually found near the top of permafrost (Smith et al., 2005;Jorgenson et al., 2006;Woo et al., 2008;Vonk et al., 2013).In these regions, thawing permafrost can result in ground ice erosion and displacement of sediments, carbon and nutrients by drainage (Rowland et al., 2010;Godin et al., 2014;Harms et al., 2014).This thermo-erosion process has especially been observed across North-America (Grosse et al., 2011), in Siberia (Günther et al., 2013) and in the Antarctic Dry Valleys (Levy et al., 2008).On Bylot Island, Nunavut, thermo-mechanical erosion by water has initiated permafrost tunneling and the development of gully networks in aeolian, organic and colluvial depositional environments of nearly 158 000 m 2 (Fortier et al., 2007;Godin and Fortier, 2012a;Godin et al., 2014;Veillette et al., 2015).A fine-scale spatio-temporal monitoring of the six largest gullies showed that their development rates ranged from 14 to 25 m yr −1 going up to 80 m yr −1 during their inception (Godin and Fortier, 2012b), leading to substantial changes in the drainage network and increasing eroded area throughout the valley (Godin et al., 2014).To date, the contribution of thermal-erosion and lateral erosion processes in permafrost feedbacks to climate has yet to be documented.
Many observational and experimental studies have highlighted shifts in tundra plant community structure and plant species productivity in response to warming temperatures (Jonsdottir et al., 2005;Hudson and Henry, 2010;Epstein et al., 2013;Naito and Cairns, 2015).In contrast, little is known about how thermo-erosion gullying affects plant community structure and plant species abundance.Yet, this information is urgently needed as vegetation plays an important role in structuring Arctic ecosystems and regulating permafrost response to climate change (Jorgenson et al., 2010;Gauthier et al., 2011;Legagneux et al., 2012).Wetlands serve as preferred grounds for Arctic herbivores such as snow geese (Gauthier et al., 1996;Massé et al., 2001;Doiron et al., 2014).They are also expected to produce more methane compared to shrub-dominated areas (Olefeldt et al., 2013;Nauta et al., 2015;Treat et al., 2015).
The present study aimed at examining plant community patterns following thermo-erosion gullying processes.Bylot Island, where geomorphological and ecological processes in response to climate change have been monitored for over 2 decades (Allard, 1996;Fortier and Allard, 2004;Gauthier et al., 2013;Godin et al., 2014), offered a unique opportunity to specifically assess the response of wetlands to gullying.The following questions were addressed: (1) to what extent thermo-erosion gullying modifies environmental conditions of low-centered wetland polygons?(2) How do plant communities cope with these geomorphological changes, i.e. do we observe shifts in plant diversity, abundance and productivity?

Study area
This study took place in the Qarlikturvik valley of Bylot Island, Nunavut, Canada (73 • 09 N, 79 • 57 W; Fig. 1a).Bound to the North and South by plateaus < 500 m a.s.l., it connects C-79 and C-93 glaciers to the Navy Board Inlet sea via a proglacial river.The sampling sites were specifically located on the valley floor (ca.65 km 2 ), which is characterized by a low-centered polygon landscape that has resulted from ice wedges development and sediment accumulation during the late Holocene (Ellis and Rochefort, 2004;Fortier and Allard, 2004;Fortier et al., 2006;Ellis et al., 2008).Two baseline vegetation types can be recognized.Wetlands, often associated with intact low-centered polygons, represent ca 23 % of the valley area (Hughes et al., 1994) and are dominated by sedges (Carex aquatilis, Eriophorum angustifolium, Eriophorum scheuchzeri), grasses (Dupontia fisheri, Pleuropogon sabinei; Gauthier et al., 1995) and fen mosses (Drepanocladus spp.; Ellis et al., 2008;Pouliot et al., 2009).Mesic environments, such as low-centered polygon rims, gently sloping terrain and hummocky tundra, support a more diverse group of species including Salix spp., Vaccinium uliginosum, Arctagrostis latifolia, Poa arctica and Luzula confusa with Aulacomnium spp. as dominant moss species (Zoltai et al., 1983).As a result, the valley houses many herbivores such as snow geese (in summer) and brown and collared lemmings (Gauthier et al., 1995;Gruyer et al., 2008), thus representing a critical environment for tundra food web (Gauthier et al., 2011;Legagneux et al., 2012).

Field sites
Our work was specifically conducted around three gullies that were selected among the 36 identified in the valley (Godin and Fortier, 2012b).These gully networks have originated from snowmelt water infiltration into cavities of the frozen active layer and the subsequent formation of underground tunnels that have ended up collapsing (Fortier et al., 2007;Godin and Fortier, 2010).The gullies R08p and R06, respectively 835 and 717 m long, are characterized by ongoing thermo-erosion processes (Fortier et al., 2007;Godin and Fortier, 2012b) whilst the gully RN08, 180 m long, has not been actively eroding in recent years.A total of 197 sampling sites were randomly selected around the three gullies (Table 1; Fig. 1b) and classified into one of four categories (referred hereafter as habitats) that represented the two baseline vegetation types (wet and mesic) as well as increasing levels of disturbance related to thermo-erosion processes.The disturbed habitats were sorted via a visual assessment of the low-centered polygon rim integrity coupled with a recent close monitoring of drainage system development along the gullies (Fortier et al., 2007;Godin and Fortier, 2012a, b).The habitats were defined as follows: (i) intact low-centered wet polygons (n = 62) that were not affected by gullying.Their elevated rims enclose a central depression that retains snow cover during winter and is flooded by snowmelt water during spring (Woo and Young, 2006;Minke et al., 2007).These polygons are hydrologically independent, with their water content representing the balance between precipitation inputs (snow and rain) and evapotranspiration outputs (Fortier et al., 2006); (ii) less than 5-year (recently) disturbed polygons (n = 44), located along the most recent sections of the gullies, with partially degraded rims, incomplete drainage and heterogeneous water content; (iii) more than 5-year (older) disturbed polygons (n = 43), with heavily breached rims adjacent to the gully channels and substantial or complete drainage; (iv) mesic environments (n = 48), with distinct heterogeneous mesic vegetation.They are found on the rim of polygons and in adjacent areas and are generally not induced by thermo-erosion gullying but rather dependent on local sedimentary and hydrological dynamics.

Environmental condition monitoring
Daily precipitation was recorded with a manual rain gauge throughout summer 2010 at the base camp, located 700 m west of the gully R08p (Gauthier et al., 2010).Soil (top 10 cm) moisture was recorded at the center of each sampling site using ECH 2 O EC-5 moisture sensors (accuracy of ±3 % VWC, resolution of ±0.1 % VWC) connected to Em5b data loggers (Decagon Devices).Measurements were carried out in 2010, both 5 July (early season) and 30 July (late season) along the gullies R08p and R06, and only 30 July along the gully RN08.Thaw-front depth was recorded at the center of each sampling site using a graduated steel rod driven into the thawed active layer.The data collection spanned 2 years with measurements in July 2009 and 2010 in the polygons of the gullies R08p and R06 and in July 2010 in those situated along the gully RN08.One measure was taken in each of the wet and disturbed sampled polygons whereas three random measurements were conducted in each mesic environment site because of the heterogeneity of this type of habitat.In this case, a mean per site was calculated prior to analyses.

Plant community characterization
Species richness and abundance were determined at each site in July 2009 or 2010 using three randomly placed 70 cm × 70 cm quadrats that were vertically photographed at ca 1.3 m from the ground (see detailed protocols in Chen et al., 2010 andthe IPY CiCAT, 2012;Fig. 2).Abundances of vascular plants, lichens, mosses, Nostoc spp., fungi, cryptogamic crust, bare ground, litter, vascular plant standing dead, standing water, signs of grubbing and goose feces were evaluated as cover percentages using photography analyses (Perreault, 2012).Daubenmire cover abundance classes (Daubenmire, 1959) were used on each quadrat picture overlain by a 7 cm grid to evaluate species cover as the projection on the ground of all species above-ground parts.
Five sampling sites per habitat were also randomly selected along the gullies R08p and R06 to measure aboveground biomass of graminoid species.At each site, an exclosure of 1 m × 1 m was made of chicken wire 30 cm high and supported by wooden stakes at each corner (see Gauthier et al., 1995).Exclosures were set up in early July 2012 to avoid any significant grazing by geese.Above-ground biomass was harvested inside the exclosures near peak production in early August 2012 using random grids of 25 cm × 25 cm for wet and recently disturbed polygons and of 50 cm × 20 cm for older disturbed polygons and mesic environments.Two different grid sizes were used because of the difference in structure of the vegetation (herbaceous vs. shrubs) associated with the habitat heterogeneity (Legagneux et al., 2012).All graminoids present in the random grids were cut to a standard height, i.e. at an average of 1 cm below the moss surface (Gauthier et al., 1995;Doiron et al., 2014), and live biomass was sorted as follows: Carex aquatilis, Eriophorum angustifolium and Eriophorum scheuchzeri (Cyperaceae); Luzula arctica and Luzula confusa (Juncaceae); Anthoxanthum arcticum, Arctagrostis latifolia, Dupontia fisheri and Festuca brachyphylla (Poaceae).Biomass was then oven-dried at 65 • C until constant dry weight and further weighed to ±0.0001 g using an electronic weighing scale.

Statistical analyses
Differences in soil moisture, thaw-front depth and graminoid above-ground biomass among habitats were tested with a generalized linear mixed model (procedure MIXED, REML method in SAS, version 9.4, SAS Institute, Cary, NC, USA).Soil moisture, thaw-front depth as well as date or year of measurements and the interaction terms were treated as fixed factors and gully as a random factor.Type III sums of squares were used for the calculation of fixed effect F statistics while random effects were assessed using a log likelihood ratio test from the full and reduced models (Littell et al., 2006).Post hoc contrasts were performed to ascertain specific differences between habitats at alpha < 0.05 using the LSMEANS statement and Bonferroni adjustment.Canonical correspondence analyses (CCA) were conducted to test unimodal relationships between habitats and environmental variables (ter Braak, 1986;Zuur et al., 2007) using the "vegan" package (Oksanen et al., 2015) in R 3.1.3(R development core team, 2015).Two matrices were elaborated: one of 197 sites × 65 taxa using their mean cover per site, and another of 197 sites × 8 environmental characteristics comprising the following continuous variables: soil moisture, thaw front depth and cover estimates of litter, bare ground, grubbing, vascular plant standing dead, standing water and goose feces.Soil moisture measurements obtained on 30 July 2010 were used in the analyses.

Environmental conditions
In 2009, above-average spring temperatures led to a rapid snowmelt (16 June) while summer was one of the driest on record (Gauthier et al., 2009).In 2010, despite a relatively warm spring (0.26 • C above normal), the high snow pack at the end of the winter (41.6 cm on 31 May) delayed the snowmelt in the lowlands to 28 June, which was a week later than normal.Summer was characterized by warm and sunny conditions as well as below-average precipitations (cumulative rain of 84 mm vs. long-time average of 92 mm; Gauthier et al., 2010).However, the 36 mm received in only 5 days in mid-July 2010 significantly increased soil moisture between the two monitoring dates in all habitats (significant date effect: df = 1, F = 88.99,P < .001;Fig. 3) .There was a significant difference overall in soil moisture among habitats (df = 3, F = 79.86,P < 0.001), which was associated with differences between wet polygons and the other habitats (5 July: df = 3, F = 33.41,P < 0.001, 30 July: df = 3, F = 47.36,P < .001;Fig. 3, Appendix Table A1).Throughout the summer, soil moisture was approximately 40 % higher in wet polygons compared to disturbed polygons and mesic environments.We also found a significant difference in thaw-front depth between wet polygons and the other habitats (2009: df = 3, F = 21.30,P < 0.001, 2010: df = 3, F = 33.86,P < 0.001; Fig. 4, Table A1).Thaw-front depth of wet polygons was approximately 10 cm deeper than in disturbed polygons and mesic environments for both years of monitoring (35-36 vs. 25-27 cm).The two types of disturbed polygons thus showed similar soil moisture and thawfront depth despite their differing time since disturbance inception.There was no significant effect of gully location (df = 2, LLR = 1.6, P = 0.21 for soil moisture, df = 2, LLR = 0.0, P = 1.0 for thaw-front depth).

Plant community characterization
A total of 18 vascular plant families encompassing 59 species were sampled throughout the study (Table A2).The greatest species richness was found in polygons that were disturbed for at least 5 years and where both hydrophilic and mesic species were present (Table 2).The transition from wet polygons to mesic environments was accompanied by significant changes in vascular plant community composition, especially with the decline in Cyperaceae and Poaceae cover and the emergence of Salicaceae species (Table 2).Carex aquatilis and Dupontia fisheri were respectively present in 100 and 93 % of the wet polygons sampled where they accounted for 52 and 26 % of the total vascular plant cover.A1 for sample sizes and post-hoc contrasts.* * * P < 0.001, ns: statistically non-significant effect.
They were found in only 47 and 16 % of mesic environments accounting for 9 and 0.45 % of the total vascular plant cover.In contrast, Salix arctica and Arctagrostis latifolia, which were found in approximately half of the wet polygons accounting for 3 % of the total vascular plant cover, were present in 98 % of the mesic environments where they respectively accounted for 50 and 14 % of the total vascular plant cover.Differences among habitats were also noted in non-vascular taxa.Abundance of lichens such as Cladonia spp., Stereocaulon spp.and Peltigera spp.increased in polygons disturbed for at least 5 years and mesic environments (Table 2).Mosses were mostly found living in wet polygons and mesic environments and dried (i.e.dead) in disturbed polygons (  Drepanocladus species replaced by Aulacomnium species in mesic environments.Moreover, we observed vegetation changes through the decline of graminoid above-ground biomass which varied significantly among habitats (df = 3, F = 11.59,P < 0.001; Fig. 5a; Table A1).Graminoid biomass was nearly five times greater in wet than in mesic environments (29.2 vs. 5.9 g m −2 ) and decreased twofold between < 5-year disturbed and > 5-year disturbed polygons (22.3 vs. 9.3 g m −2 ; Figs. 2 and 5a).Differences were mainly driven by the decline of hydrophilic species, i.e.Carex aquatilis, Eriophorum scheuchzeri, Anthoxanthum arcticum and Dupontia fisheri, between wet and mesic habitats (28.8, 19.7, 3.6 and 2.5 g m −2 in wet, < 5-year disturbed, > 5-year disturbed and mesic polygons, respectively; Fig. 5b).Above-ground biomass of Luzula spp., Arctagrostis latifolia and Festuca brachyphylla was contrastingly nine times greater in mesic than in wet habitats (3.30 vs. 0.35 g m −2 ) and twenty-four times greater in > 5-year disturbed than in < 5-year disturbed polygons (5.61 vs. 0.23 g m −2 ).

Relationships between plant communities and environmental variables
The first two axes of the canonical correspondence analysis retained 14 % of the vegetation data variance and 80 % of the vegetation-environment relationship variance (Table 3).Five of the eight environmental variables tested were significant within the canonical model (P < 0.05, 999 permutations), but only three -litter cover, thaw-front depth and soil moistureshowed high correlations with the canonical axes (Table 4).Thaw-front depth and soil moisture were strongly related to the first axis, while litter cover was mainly associated with the second axis (Table 4).Altogether, these variables dis- criminated well the four studied habitats of the Qarlikturvik valley.Wet polygons were mainly related to high soil moisture and substantial thaw-front depth whilst mesic environments were associated with greater litter cover (Fig. 6).The gradual vegetation transition was also observed from < 5year to > 5-year disturbed polygons following the soil moisture shift in these habitats (Fig. 6).

Discussion
Sustainability of wetlands at high latitudes essentially relies on perennial frozen ground that prevents drainage and allows wet soil conditions (Woo and Young, 2006;Ellis et al., 2008).However, snowmelt water run-off through ice-wedge polygon landscapes can initiate thermal erosion of the permafrost and the development of gullies (Fortier et al., 2007;Godin and Fortier, 2014).We showed here that permafrost gullying significantly altered wetlands by changing the orig-  inal polygon microtopography, and decreasing soil moisture and thaw-front depth of disturbed polygons along the gullies.Vegetation was sensitive to this process, and mesic environment plant species gradually replaced hydrophilic species within 5 to 10 years, although the full transition has yet to be reached.This vegetation turn-over can have substantial consequences on wildlife biology, permafrost stabilization and ecosystem-level greenhouse gas emissions (Blok et al., 2010;Doiron et al., 2014;M. T. Jorgenson et al., 2015;McEwing et al., 2015).

Transient environmental conditions
Thermo-erosion gullying led to a significant decrease in soil moisture and thaw-front depth of breached polygons.Both older and recently disturbed polygons had similar soil moisture and thaw-front depth while differing in time since disturbance, which shows that the change in polygon environmental conditions after permafrost disturbance was rapid.The decrease in soil moisture following polygon rim erosion is consistent with what has been previously observed in gullied areas (Seppälä, 1997;Godin and Fortier, 2012a, 2014, 2015;Harms et al., 2014) and concurs with a modeling analysis showing that the transformation of low-centered to highcentered polygon landscape following ice wedge degradation is accompanied by a significant alteration in the water balance partitioning (Liljedahl et al., 2012).In our study, all types of polygons were recharged by snowfall and summer rainfall, yet disturbed habitats had lower soil moisture than wet polygons and a thorough examination of moisture evolution throughout an entire summer showed that soil moisture of breached polygons was significantly more variable than that of wet polygons at both intra-and inter-polygonal scales (Godin et al., 2015).Given that soil moisture is an important driver of plant community composition (Dagg and Lafleur, 2011), it is no surprise that we observed a shift in vegetation following changes in moisture regime.Decreasing soil moisture in the center of disturbed polygons came with decreasing thaw-front depth, which was expected given that active layer thickness is closely related to soil moisture (Nelson et al., 1999;Hinzman et al., 2005;Minke et al., 2009;Wright et al., 2009;Gangodagamage 4).Standing dead represents the cover of dead attached vascular plants.et al., 2014).This result, however, contrasts with the active layer thickening generally observed in response to climate warming (Tarnocai et al., 2004;Woo et al., 2007;Akerman and Johansson 2008;Smith et al., 2009;Nauta et al., 2015), and this is likely due in part to ground surface subsidence and drainage which follows ice-rich permafrost thawing (Shiklomanov et al., 2013) and in part to snow accumulation patterns (Godin et al., 2015).Within 5 years of drainage, thaw-front depth in disturbed polygons decreased by 37 % compared to that in wet polygons.This is mainly explained by heat capacity of water and the higher thermal conduction rates in wetter polygons that provide substantial melt energy to the frost table (Nelson et al., 1997;Hinzman et al., 2005;Wright et al., 2009;Romanovsky et al., 2010).This effect is also sharpened by the low thermal conductivity of drier moss carpets in disturbed habitats (Wright et al., 2009) and reduced local snow conditions within the polygons adjacent to the gullies (Godin et al., 2015).

Vegetation changes
Overall, the floristic composition of our sampling sites is in line with previous field surveys conducted in the same area (Gauthier et al., 1996;Duclos, 2002;Doiron, 2014).The presence of Carex aquatilis, Eriophorum scheuchzeri and Dupontia fisheri characterizes well the typical vegeta- tion of Arctic wetlands (Jorgenson et al., 2013;Sandvik and Odland, 2014;Lara et al., 2015) whilst that of Arctagrostis latifolia, Luzula and Salix spp.are common features of Arctic mesic environments (Audet et al., 2007;Sjögersten et al., 2008).Disturbed polygons were the most diverse habitats given that they offered a middle-range state between wet and mesic conditions where hydrophilic species were still present while mesic environment ones had successfully established.The development of gullies in the Qarlikturvik valley and the subsequent drainage of adjacent low-centered polygons have led within 5 to 10 years to a gradual change in plant communities with vegetation of disturbed polygons leaning toward a new equilibrium, that of mesic environments.Mesic environment species such as Luzula and Salix spp.have established or increased in cover following the decrease in soil moisture and thaw-front depth and replaced hydrophilic Cyperaceae and Poaceae.The secondary succession pioneered here by the gullying process in disturbed polygons follows the directional-species replacement model examined by Svoboda and Henry (1987).However, by occurring within 5 to 10 years, it has been remarkably more rapid than what is usually documented for the High Arctic where perennial plant communities are largely resistant to disturbance (Hollister et al., 2005;Jonsdottir et al., 2005;Hudson and Henry, 2010) and succession dynamics are slow due to short growing seasons and low summer temperatures (Svoboda and Henry, 1987).For instance, plant cover of northeastern Alaska changed little over a 25-year period despite a significant rise in summer temperatures (J.C. Jorgenson et al., 2015).This gradual yet rapid species replacement has been triggered in our system by the hydrological and thermal shift caused by gullying and favored by the mosaic of wet and mesic habitats allowing for a substantial pool of species with both vegetative and sexual reproduction.
In the canonical ordination analysis, the soil moisture gradient discriminated wet polygons from the other habitats as well as recently disturbed from older disturbed habitats.The 37 % decrease in soil moisture between wet and disturbed polygons represents a drastic change of conditions for plant communities and is of similar magnitude than what has been documented in Alaskan drying wetlands as a result of increasing temperatures (Klein et al., 2005).The strong influence of soil moisture in separating plant community types at high latitudes has indeed been well documented (Hinzman et al., 2005;Daniëls and de Molenaar, 2011;Daniëls et al., 2011;Sandvik and Odland, 2014).Four other variables significantly influenced the distinction among habitats: (i) thaw-front depth discriminated habitats in the same direction than soil moisture with a 30 % decrease in disturbed polygons and mesic environments compared to wet polygons, which was expected since these two factors are closely related (see Sect. 4.1); (ii) litter cover separated mesic environments from the others, which may be explained by increased organic matter related to greater shrub abundance in mesic environments (Zamin et al., 2014); (iii) vascular plant standing dead separated wet and recently disturbed polygons from the other habitats, which can be explained by the senescence of Cyperaceae tillers that are highly abundant at these locations (Fig. 5); (iv) goose feces were mainly associated with older disturbed and mesic environments.While this may suggest a higher use of these habitats by geese, the slower degradation of feces in dryer habitats cannot be ruled out; this has yet to be tested.
The shift in vegetation composition in disturbed polygons was accompanied by significant changes in biomass.Aboveground biomass of graminoids was the greatest in wet polygons, which is concordant with the fact that wetlands are the most productive habitats of forage plants in the Arctic (Sheard and Geale, 1983;Duclos, 2002;Doiron, 2014).It gradually decreased in disturbed polygons as conditions became closer to those of mesic environments.Compared to the immediate change in environmental conditions, we nonetheless observed a lag-time in vegetation response to thermo-erosion related disturbances as graminoid biomass differed significantly between recently and older disturbed polygons.In our study, graminoid above-ground biomass of wet polygons was 35 % lower than what Cadieux et al. (2008) found via a long-term plant monitoring on By-lot Island (45.2 g m −2 ), and 62.7 % lower than what Gauthier et al. (2012) measured in the most productive wetlands of the Qarlikturvik Valley (78.4 ± 10.5 g m −2 ).These contrasts may be explained by varying species composition and to a lesser extent by earlier plant harvesting in our case.Indeed, while we focused on wet polygons dominated by Carex aquatilis, Cadieux et al. (2008) and Gauthier et al. (2012) worked on wet polygons dominated by Dupontia fisheri and Eriophorum scheuchzeri.Because our study was part of a large-scale multisite project on wetland carrying capacity for snow geese (Legagneux et al., 2012;Doiron, 2014), we only focused on forage plant (i.e.graminoids) biomass and did not sample forbs or shrubs.Since above-ground biomass of graminoids account for more than 90 per cent of vascular plant biomass in wetlands (Gauthier et al., 1995), we provide here an accurate estimate of the total above-ground biomass that can be found in these habitats.However, the total aboveground biomass in mesic environments was probably underestimated.For instance, biomass of shrubs and forbs respectively ranged between 22 and 48 g m −2 and between 6 and 20 g m −2 in the mesic environments adjacent to our study area (E.Lévesque, unpublished data).Overall, total aboveground biomass in wetlands and mesic environments is of similar order of magnitude (50.5 g m −2 ± 2.8 SE in wetlands and 44.2 g m −2 ± 6.8 SE in mesic tundra for the period 2007-2009; Legagneux et al., 2012).

Impacts on ecosystems
It is likely that the replacement of hydrophilic plants by mesic vegetation will severely impact wildlife biology.The Qarlikturvik valley of Bylot Island represents an important summer habitat for greater snow geese (Legagneux et al., 2012).It is well documented that this species mostly relies on wetlands for food resources (Gauthier et al., 1995(Gauthier et al., , 2011)), especially because graminoids are easily digested due to their low fiber concentration and rich nutritive elements (Sedinger and Raveling, 1989;Manseau and Gauthier, 1993;Audet et al., 2007).For instance, geese removed respectively 40 and 31 % of the total annual production of Dupontia fisheri and Eriophorum scheuchzeri during the period 1990-2007 (Cadieux et al., 2008).It remains to quantify the extent to which gullying alters wetland carrying capacity.In addition, the presence of ponds in wetlands provides geese refuges from predators such as the arctic fox (Hughes et al., 1994;Lecomte et al., 2009), and their disappearance might also change predator-prey interactions.
Effects of gullying-induced vegetation changes may finally be visible on variations of greenhouse gas emissions.There is evidence for a strong vegetation control on methane emission from wetlands (Olefeldt et al., 2013;Lara et al., 2015;McEwing et al., 2015;Tveit et al., 2015).In wet polygonal tundra of Northern Siberia, Kutzbach et al. (2004) found for instance that dense Carex aquatilis stands emitted more methane than sites with low Carex densities.Overall, wet-land and lake expansion are thought to increase methane emission but also carbon storage (Myers-Smith, 2005;Nauta et al., 2015;Treat et al., 2015;Bouchard et al., 2015).We can therefore expect that the reverse transition from wet to mesic environments observed within our low-centered polygon landscape would lead to reduced methane emission and increased CO 2 emission through enhanced decomposition.However, no general pattern on ecosystem responses to decreased water table position and subsequent gas emissions has emerged to date (see Grosse et al., 2011 for review).It will thus be crucial to determine in the near future the specific evolution of Salix and Luzula spp.primary production in mesic environments in order to accurately predict the effects of wetland retreat on methane and soil organic carbon cycles.

Conclusions
This study illustrates that changes in the hydrological and thermal regimes following thermo-erosion gullying processes boost landscape transformation from wet to mesic habitats, providing evidence that permafrost disturbance is a critical component of ecosystem modification at high latitudes.Ecological studies should consequently start using an approach that integrates disturbed permafrost monitoring if one wants to more efficiently document climate change effects on arctic terrestrial ecosystems.In addition, our latest field observations showed that hydrology and thaw regimes of breached polygons have yet to reach equilibrium with new conditions.Similarly, vegetation remains in transition given that, 10 years after disturbance, the cover of dominant shrubs and mesic bryophytes in disturbed polygons is still lower than in adjacent mesic environments.It is currently not possible to predict how long these species would take to out compete declining species and cryptogamic crust and reach a new mesic environment equilibrium.This current state underscores the importance of long-term monitoring of permafrost and its associated vegetation.In addition, more work should be devoted to the feedback effects of plant communities and vegetation succession on thermal and mechanical stabilization dynamics of disturbed permafrost terrains.This is especially needed since plant community differences between disturbed and intact sites can last several centuries (Cray and Pollard, 2015).

Figure 2 .
Figure 2. The four habitat types studied in the Qarlikturvik valley, Bylot Island, Nunavut.The close view at the bottom right of each picture represents the 70 cm × 70 cm quadrats that were used to determine species richness and abundance in each sampling site.

Figure 3 .
Figure 3. Soil moisture monitored early and late July 2010 in the four habitat types studied in the Qarlikturvik valley, Bylot Island, Nunavut.The 10th percentile, lower quartile, median (solid line), mean (dashed line), upper quartile and 90th percentile are shown.See TableA1for sample sizes and post-hoc contrasts.* * * P < 0.001, ns: statistically non-significant effect.

Figure 4 .
Figure 4. Thaw-front depth monitored in July 2009 and 2010 in the four habitat types studied in the Qarlikturvik valley, Bylot Island, Nunavut.The 10th percentile, lower quartile, median (solid line), mean (dashed line), upper quartile and 90th percentile are shown.See TableA1for sample sizes and post-hoc contrasts.* * * P < 0.001, ns: statistically non-significant effect.

Table 1 .
Distribution of the sampling sites per habitat type and per gully.

Table 2 .
Species richness, family total cover and species mean cover of vascular taxa as well as mean cover of non-vascular taxa in each of the four habitat types studied in the Qarlikturvik valley, Bylot Island, Nunavut.Mean species richness is given for sampled areas of 49 dm 2 (70 × 70 cm quadrats).Numbers in brackets denote the number of species inventoried in each family.= cover < 0.01 %; < = cover < 0.1 %.Species names were retrieved (2011) from the Integrated Taxonomic Information System (ITIS) (http://www.itis.gov).

Table 3 .
Canonical correspondence analysis of the vegetation and environmental data gathered in four habitat types in the Qarlikturvik valley, Bylot Island, Nunavut.CCA-1: first canonical axis; CCA-2: second canonical axis.

Table 4 .
Canonical correspondence analysis of the vegetation sampled along three gullies in the Qarlikturvik valley, Bylot Island, Nunavut, and biplot scores for environmental variables.CCA-1: first canonical axis; CCA-2: second canonical axis.Statistically significant values (P < 0.05) after 999 permutations are shown in bold.Standing dead represents the cover of dead attached vascular plants.