Importance of Marine-Derived Nutrients Supplied by Planktivorous Seabirds to High Arctic Tundra Plant Communities

We studied the relative importance of several environmental factors for tundra plant communities in five locations across Svalbard (High Arctic) that differed in geographical location, oceanographic and climatic influence, and soil characteristics. The amount of marine-derived nitrogen in the soil supplied by seabirds was locally the most important of the studied environmental factors influencing the tundra plant community. We found a strong positive correlation between δ15N isotopic values and total N content in the soil, confirming the fundamental role of marine-derived matter to the generally nutrient-poor Arctic tundra ecosystem. We also recorded a strong correlation between the δ15N values of soil and of the tissues of vascular plants and mosses, but not of lichens. The relationship between soil δ15N values and vascular plant cover was linear. In the case of mosses, the percentage ground cover reached maximum around a soil δ 15N value of 8‰, as did plant community diversity. This soil δ15N value clearly separated the occurrence of plants with low nitrogen tolerance (e.g. Salix polaris) from those predominating on high N content soils (e.g. Cerastium arcticum, Poa alpina). Large colonies of planktivorous little auks have a great influence on Arctic tundra vegetation, either through enhancing plant abundance or in shaping plant community composition at a local scale.


Introduction
Polar terrestrial ecosystems are distinctive, as they develop and function under very harsh conditions [1]. Among the main recognized ecological factors affecting the development and dynamics of terrestrial vegetation in polar regions are: air temperature [2,3], soil moisture [4,5], soil pH [5], nutrient availability [6], snow cover [7], proglacial chronosequences [8], dispersal limitations [9] and natural disturbances [4]. These factors are modulated at the macroclimatic scale, e.g. resulting from geographical separation, as well as by microclimatic features such as topography (elevation and exposure) [10] and, in turn, influence the structure and distribution of polar vegetation [11,12].
The tundra ecosystem is generally poor in nutrients [1]. However, at a local scale, marine birds and mammals deposit large amounts of excrement (rich in nitrogen, phosphorus and other elements) near their breeding colonies or permanent resting sites, fertilizing these areas [13][14][15]. Nutrients delivered in this way play an important role in defining local soil chemical characteristics, enhancing vegetation productivity and biomass, and influencing the spatial distribution of plants [13][14][15][16][17]. Ornithogenic compounds assimilated by plants are subsequently transferred to higher trophic levels and finally, through decomposition, back to soil [18,19].
Little auks (dovekies, Alle alle) play a potentially crucial role in Arctic tundra development. They are the most numerous and widespread High Arctic seabird, with a global population of ca. 37 million pairs [20]. The Svalbard population alone has been estimated to comprise more than 1 million pairs [20]. Their role as biovectors is that they forage in the sea and deposit droppings on land. The food of little auks is comprised mostly of copepods (Calanus glacialis and C. hyperboreus in cold Arctic waters, and C. finmarchicus in warmer Atlantic waters) [21]. During the breeding season, little auks supply up to ca. 60 t km -2 dry mass of faecal matter in the vicinity of their colony in Hornsund (south-west Spitsbergen) [22]. Our previous study [15] showed that such large levels of nutrient deposition have strong influence on soil physical and chemical parameters there, and also create a steep fertility gradient from the colony to the sea.
It is well known [23,24] that nitrogen and carbon play fundamental roles in the growth and development of plants, and are amongst the most important resources in biogeochemical cycles. Depending on whether these elements are of marine or terrestrial origin, the ratios of their stable isotopes ( 15 N/ 14 N expressed as δ 15 N, and 13 C/ 12 C as δ 13 C, respectively) are different [23]. Since the proportion of heavier isotopes is generally greater in marine organic matter, and further increases when passing through successive trophic levels, in this study we used δ 15 N and δ 13 C signatures of soil and plant tissues as indicators to identify seabird influence [25,26].
Although the effects of seabird colonies on tundra vegetation have been widely reported [16,18,27,28], there are no quantitative studies comparing the combined response (competitive ability) of specific plant species, or the entire tundra plant community, to multiple environmental factors affecting them at the same time. In this study we tested the relative importance of the following factors for tundra development: (1) total nitrogen and carbon content and isotopic signatures of soil, (2) soil characteristics (moisture, pH, conductivity), and (3) geographical location.
We test the hypothesis that nutrients derived from the marine ecosystem via colonial seabirds are locally the main factor determining the abundance and structure of plant communities in Svalbard. The other aims of this study were also to identify plant species most dependent on birds' nutrient enrichment and to describe their response along the fertilization gradient.

Study area
The study was conducted at five locations within the Svalbard archipelago, differing in local climate regimes and the sizes of local little auk populations (Table 1, Fig 1).

Data sampling
The study was conducted in July and August, during expeditions to Hornsund (2005Hornsund ( , 2007, Magdalenefjorden and Aasefjelet (2007 and 2009), Bjørnøya (2008), and Isfjorden (2010). Samples were collected along 12 line transects established in the five locations described above. The upper parts of the transects consisted of vegetation-covered (to differing extents) rock debris, while the lower parts, approaching the sea/lake shore, were flat and more waterlogged. Each transect, depending on its local geomorphology condition, consisted of 5-10 plots (160×160 cm each) that were located from the respective transect's starting point (plot 1) as follows: plot 2 (6 m), 3 (15 m), 4 (29 m), 5 (49 m), 6 (79 m), 7 (125 m), 8 (193 m), 9 (296 m), 10 (449 m), 11 (680 m), and 12 (1026 m), as described by Zwolicki et al. [15]. Within each sampling plot we assessed: (i) vegetation percentage cover, and collected (ii) three soil cores from three sites along the same diagonal of each plot (from the centre and from the two corners of each square), for physical, chemical and isotopic analyses, and (iii) three samples of taxa representing monocots, eudicots, mosses and lichens, for isotopic analyses.

Guano deposition measurements
Guano deposition measurements were obtained in summer 2005 in Hornsund location only, during the study which compared the influence of guano deposition on soil chemistry between plankton-eater and fish-eater seabird colonies [15]. In the present study we used these data in order to reveal the relation between guano deposition rate and values of δ 15 N in the soil.
Guano deposition was assessed using black plastic sheets (150×150 cm) placed next to each sampling plot, along the colony and control transects. Exposure time depended on weather conditions, with the range of 20-36 hours. After exposure, a photograph of each sheet was taken and the sheets were cleaned for re-exposure. The area covered with birds' droppings in each photograph was measured using SigmaScan Pro 5.0.0 software. Situated close to the Atlantic water masses, but also partly influenced by cold Arctic waters; represents typical maritime climatic conditions.
In order to provide repeatable estimates of guano deposition (dry mass) from the photographs, we performed an initial calibration by exposing stiff plastic sheets (150×150 cm), covered with very thin plastic film of a known mass. After these sheets were exposed and photographed, the plastic films were collected, dried and re-weighed to obtain the dry mass of droppings and to calculate a regression equation (y = 0.003 x, R 2 = 0.7, n = 31; where: x-area covered by guano (cm 2 ), and y-guano dry mass (g)) [15].

Vegetation abundance and species composition
Within each sampling plot (total n = 97) we identified vascular plant species, and visually estimated the percentage contributions of individual species, and also the contributions of vascular plant and moss cover to the total vegetation cover. In order to describe changes in species αdiversity along the ornithogenically-influenced gradients, we decided to calculate Hill's N2 diversity index, which measures the effective number of species and is linearly related to the number of species [29], separately for each plot (S1 File).

Physicochemical soil properties
Each sample was collected with a shovel from the soil surface layer (to a depth of ca. 5-10 cm) and contained about 500 cm 3 of soil. At sampling sites with very compact vegetation, we removed and discarded the upper layer of live and dead, poorly decomposed, plant material. Immediately after collection of samples (maximum 24 hours), in the field laboratory or in University of Svalbard, we took two sub-samples, 80 cm 3 each, from each sample and performed assessment of: 1. soil moisture (%)-Soil sub-samples were weighed with electronic scales (precision 0.1 g) before and after oven-drying (60°C) to a constant mass. Soil moisture was defined as: soil moisture = ((wet mass-dry mass) dry mass -1 ) Ã 100%.
2. soil conductivity (μS cm -1 ) and pH-Soil samples of 80 cm 3 were mixed with 160 cm 3 of distilled water. The mixture was shaken for ca. 20 min and then filtered through a sieve (0.5 mm mesh). Conductivity and pH were quantified in the filtrate using a pH/conductivity/ salinity meter CPC-401 (Elmetron) (S2 File).

Stable isotope analyses
To assess δ 15 N and δ 13 C signatures in soil the remaining part of each soil sample was sieved (0.25 mm mesh) to remove stones and large plant debris, and ground with a vibrating mill (LMW-S, Testchem) to a grain size of less than 0.03 mm.
In case of plant and lichen tissues, we collected three samples from the aboveground parts of common vascular plants, mosses and lichens from each sample plot (not less than 5 mg dry mass in each sample). Immediately after collection they were manually cleaned of contaminants such as guano, soil particles etc., dried at 40-60°C to a constant mass and ground with a vibrating mill.
Prior to isotopic analyses we removed inorganic carbon (by adding 1 mL HCl 1N per 100 mg of soil) and lipids (using 4 ml of cyclohexane per 50 mg of soil). After this, a small amount of each sub-sample (1-2 mg, weighed with a microbalance, precision 0.001 mg) was packed into a tin capsule.
Nitrogen and carbon isotope ratios were determined by a continuous flow mass spectrometer (Thermo Fisher, Delta V Advantage) coupled to an elemental analyser (Thermo Fisher, Flash EA 1112). All samples for stable isotopes have been analyzed in one laboratory at the University of La Rochelle (France). Results were expressed in the conventional δ 15 N and δ 13 C notation, according to the equation: δ X = (R sample R standard -1 -1) 1,000 (‰), where R sample was the stable isotope ratio 15 N/ 14 N or 13 C/ 12 C in the analysed sample, and R standard was the stable isotope ratio 15 N/ 14 N or 13 C/ 12 C (respectively) in the reference material i.e. atmospheric N 2 for nitrogen and PeeDee belemnite for carbon [23].
Within all 97 sample plots from five locations, we collected 1701 soil and plant samples for stable isotope analysis. Total nitrogen and carbon content (%) was also measured in each soil and vegetation sample during the isotopic analyses (S3 File).

Statistical analyses and data management
In order to test the relationships between guano deposition and soil N and C contents and isotopic signatures, the non-parametric Spearman's rank correlation was used (due to non-normal distributions of data and a relatively low number of sampling plots per group tested). This analysis was performed only in Hornsund location to validate that the value of δ 15 N was a good predictor of the ornithogenic fertilization (S4 File).
To explore individual soil and vegetation parameters and species response to soil δ 15 N we employed General Linear Models (GLM) to illustrate simple linear relation or Generalized Additive Models (GAM) if they better described species response curves than GLM. To find the best fit, Akaike Information Criterion (AIC) was performed, using Canoco 5.03 [30]. Models with scatter plots are available in supplementary materials (S1 Fig).
To explore the influence of theoretical environmental gradients in the data, and for comparison with constrained models (with the variability explained by specific environmental factors), unconstrained models were used (PCA-Principal Component Analysis-for N and C contents, and isotopic signatures of vascular plant, moss and lichen tissues; DCA-Detrended Correspondence Analysis-for plant community composition).
Depending on the length of the gradient obtained in DCA, measured in standard deviation (SD) units, Canonical Correspondence Analysis (CCA; when SD >3) or Redundancy Analysis (RDA; SD <3) were used to examine the influence of soil variables and the geographic location factor on vegetation properties [30]. After CCA or RDA, a Monte Carlo test with 499 permutations was performed to identify which of the factors significantly influenced the model. The efficiency of the environmental factor(s) in explaining the non-random variability existing in the data (%) was calculated by dividing the percentage variability explained by a given environmental factor by that explained by the first four axes in a PCA or DCA [30].
To calculate the unique contribution of seabird-derived matter in explaining the variability of vascular plant community composition as distinct from the effects of geographic location we used Variation Partitioning [31] based on two unimodal models: (1) with the four soil properties found to be significant in the previous constrained model, i.e. soil δ 15 N, moisture and pH, and plot order, and (2) with the five geographical locations. Four degrees of freedom (df) were used in each group to stabilize the model when the Monte Carlo permutation test was performed. For multiple comparisons we used Holm's correction to control the family-wise type I errors [32].
To explore significant positive and negative relationships between individual plant species and soil δ 15 N values, a t-value biplot (with Van Dobben circles), which approximates the t-values of the regression coefficients of a weighted multiple regression, was generated [31].
The results were processed using the STATISTICA 9.1 package for correlations [33], and the CANOCO 5.03 package for ordination methods and regression models [30].

Guano deposition, soil N and C contents and isotopic signatures
With increasing deposition of guano in Hornsund, the percentage of total soil nitrogen and carbon increased significantly (r s = 0.70, p < 0.001; r s = 0.71, p < 0.001; respectively). In this location the positive correlation between the deposition of bird faeces and soil nitrogen stable isotope ratio (δ 15 N) was also highly significant (r s = 0.69, n = 20, p < 0.001), while it was negative for the carbon signature (δ 13 C; r s = -0.70, n = 20, p < 0.001).
In all five locations studied, total soil N significantly increased with the higher soil δ 15 N values, starting to grow rapidly from 8‰ onwards (GAM, p < 0.001, Table 2, Fig 2A).

Relationships between soil δ 15 N and vascular plant, moss and lichen N and C contents and isotopic signatures
Soil nitrogen isotope value was linearly related to δ 15 N of vascular plants and mosses (GLM, p < 0.001), with these two groups showing almost identical trends. In the case of lichens, the relationship was not statistically significant (GLM, p = 0.059; Fig 2B, Table 2). Redundancy analysis (RDA) revealed that soil δ 15 N was the only factor contributing significantly to explaining the variability of total N and C contents and isotopic signatures of vascular plants, mosses and lichens, and accounted for 30.8% of the total variation of these parameters (Fig 3, Table 3).  Table 2). A similar unimodal, significant relationship was found between soil δ 15 N content and the Hill's diversity coefficient (N2) (GAM, p = 0.013, Fig 2A, Table 2).

Responses of plant communities and individual species to soil δ 15 N values
The first four axes of DCA explained 30.9% of the total variability in vascular plant community composition. Of the eight tested soil properties, four significantly contributed to this variability (CCA; Table 3, Fig 4). The most important factor impacting community composition was soil δ 15 N, which explained 34.1% of the variability that was possible to explain (efficiency), followed by soil moisture (18.2%), plot order (9.3%), and pH (6.6%). The other four tested variables (total soil nitrogen and carbon, soil δ 13 C, conductivity) did not significantly influence community composition.
Based on the t-value biplot analysis performed for all 36 taxa of vascular plants, nine species were significantly correlated with soil δ 15 N (Fig 4). Five of them (Ranunculus pygmaeus, Oxyria digyna, Poa alpina, Cerastium arcticum, Cochlearia groenlandica) responded positively with increases in their percent cover. Negative relationships were found for four species (Saxifraga oppositifolia, Festuca rubra, Salix polaris and Equisetum boreale). These responses were particularly clear in the common species, C. arcticum, P. alpina and S. polaris (Figs 2D and 4A, Table 2).

The relative importance of soil properties and geographical location for plant community composition
Four statistically important factors influencing vascular plant community composition, i.e. δ 15 N content, soil moisture, pH and plot order (Table 3), and five geographical locations, explained in total 92.3% of the variability of the plant community. The variation partitioning test (p < 0.01) revealed that the soil properties explained 52.5%, while the location factor accounted for 63.2%, of this variability. The two groups of factors shared 23.3% of the explained variation ( Fig 5A). Soil δ 15 N was responsible for 32.4%, moisture accounted for 17.2%, pH for 7.8% and plot order for 6.5% of the variation contributed by soil properties (Fig 5B).

Discussion
One of the most important general features of polar terrestrial ecosystems is their chronically low nutrient availability to plants, through the poor quality of the soil which also spends much of the year frozen [34]. Therefore, any additional nutrient sources such as seabird guano are of Table 2. GLM and GAM models results of the responses of the tested variables to soil δ 15 N (response curves given in Fig 2).

Response variable
Model great importance for ecosystem functioning [18,28,35]. This study presents direct evidence of a strong relation between the proportion of nitrogen stable isotopes in soil and the level of seabird guano deposition, and confirms previous studies proposing that δ 15 N signature is a reliable measure of maritime nitrogen supply by seabirds [36][37][38][39][40].
The highest value of soil δ 15 N found in our study was 13.7‰, similar to that reported by Skrzypek et al. 2015 [40]. The value is considerably higher than that of the little auk faeces itself reported by those authors (8.1‰), probably as a consequence of microbial activity in the soil followed by ammonia volatilization [41]. Soil δ 15 N values exceeding 8‰ could also be influenced by other nitrogen sources present in the area, such as little auk carcasses, in which the N isotopic ratio in muscle tissue may reach 11.2‰ [42], and in feathers 11.4‰ [43].  The δ 15 N ratio was strongly correlated with total nitrogen amount in the soil, which confirms a general assumption that seabirds as a biovectors are locally the most important source of nutrients in Arctic tundra ecosystems [15]. Furthermore, the levels of δ 15 N in the tissues of vascular plants and mosses were linearly related with δ 15 N of the underlying soil. Finally, nitrogen isotopes in soil was the most important variable explained by the total amounts of N and their stable isotope signatures in vegetation. These relationships indicate that plants easily incorporate and rely on this source of nitrogen [44,45]. The slight differences between the relationships shown by vascular plants and mosses may be connected with the different N-loading strategies used by the two groups, i.e. absorption through the roots by vascular plants, and via the surface of leaves by mosses [46]. The relation was insignificant for lichens probably due to the utilization of atmospheric nitrogen fixation [47].
Our results indicate that the carbon stable isotope is not the best indicator of the influence of seabirds. We have not found any significant relations between δ 13 C values in the soil and in the plant tissues. This was also an insignificant variable in explaining plant community composition. Primarily, we expected that the content of 13 C would be greater in the matter originating from the sea [26,48]. However, many authors faced similar difficulties with interpretation of carbon isotope results (e.g. [48][49][50]).  Table 3. Total and explainable conditional effects of the environmental variables on vegetation total N and C content, δ 15 N and δ 13 C (RDA), and on vascular plant community composition (CCA).  The data obtained in this study demonstrate that tundra communities vary mostly between different geographical locations. However, at a local scale we documented that nutrient input, which little auks were primarily responsible for, created the most important gradients of all the tested variables, and explained the largest proportion of variation in the plant community. Our variation partitioning model characterized almost maximal efficiency (92.3%) which implies that it includes most of the environmental factors that influence the variability of studied tundra vegetation. The fraction shared between geographical location and soil variables probably results from the differences in colony size and fertilization level between locations.
Besides having isotopically enriched tissues, vascular plants and mosses clearly responded to the ornithogenic nitrogen input through increases in abundance. Total cover of vascular plants increased linearly with soil δ 15 N by up to 60%, when the nitrogen isotopic ratio achieved its maximum (almost 14‰). Total moss percentage cover also increased with soil δ 15 N level but reached the highest values (ca. 55%) around δ 15 N = 8‰, then began to decrease. This domination of vascular plants in highly fertilized conditions could be the result of their much faster nitrogen uptake allowing them to outcompete mosses [51]. The peak of the plant diversity (Hill's index) for the community was also very close to δ 15 N = 8‰. It is very suggestive that this tipping point is associated with the δ 15 N signature of bird faeces, determined in a previous study to be 8.1±0.5‰ in the Hornsund area [40].
This critical N isotopic level appeared to have a strong influence on the proportions of particular plant species in the community. Along the δ 15 N gradient in the soil, plant species distributions demonstrate clear niche segregation. Changes in nutrient availability should impact on interspecific competition, as species with an ability to respond rapidly to increased nutrients will have a competitive advantage over those adapted to low nutrient levels [51,52]. This should result in selection for plants with high growth rates and rapid tissue turnover but low nutrient use efficiencies (e.g. C. groenlandica, P. alpina, C. arcticum, and O. digyna in this study), and against slow-growing plants with long-lived leaves and high nutrient use efficiency (S. oppositifolia and S. polaris here) [53,54]. Field observations and manipulation experiments in the Arctic Aleutian Islands showed that some graminoid species spread rapidly in response to increasing amounts of seabird-derived nutrients, outcompeting slower-growing dwarf shrubs and forbs, and forming lush stands of grasses and sedges with reduced overall plant diversity [55]. The responses of nine vascular plant species in the current study (25% of the vascular plant taxa recorded) were significantly related to ornithogenic fertilization. S. polaris and S. oppositifolia are typical of species growing in poor nutrient conditions, with their highest percentage covers observed under the lowest δ 15 N values (from 2 to 4‰). These perennial dwarf woody shrubs are widespread in Arctic deserts, and are well adapted to unpredictable resource availability and limitation of growth [53]. They exemplify stress-tolerant plants which are able to store their restricted resources [56]. With increasing nutrient input we observed that the coverage of these plants declined, and there was a rapid increase and domination of other species, especially P. alpina (at levels from 6 to 10‰) and C. arcticum (from 8 to 13‰). P. alpina was first described as a dominant species in nitrophilous communities by Summerhayes and Elton [57], and C. arcticum is also known as a very abundant species occurring in the vicinity of bird colonies [58,59].
Many of the vascular plants that typify ornithogenic tundra [59] are reported to be important food for large herbivores (reindeer, geese, ptarmigan, on Svalbard). These fast-growing, nutritious and palatable plant species create large grazing areas of good quality in the vicinity of seabird colonies [60,61]. Previous studies have also shown that areas without seabirds are characterized by low-productivity tundra dominated by dwarf shrubs and forbs [18,55,62]. Similarly, on sub-Antarctic islands, nutrient subsidies from seabirds encourage a change from carpets of ferns and mosses to lush tussock grasslands [37,63]. Such selection for fast-growing graminoids in response to high nitrogen availability is typical for cold regions where annual productivity is low [26,64]. In our study, beside the grass species, the ornithogenic N-input favours species such as C. groenlandica and O. digyna, that probably do not rely exclusively on fungal or bacterial symbionts for nutrient acquisition [65,66]. As with the graminoids, both these species are characterized by high nutrient contents and relative growth rates, as well as rapid maturation and high fecundity [14].
Nutrient enrichment derived from colonial seabirds typically enhances above-ground biomass, while also locally reducing plant diversity if fertilization is considerable [27,67]. We observed this phenomenon in the areas of the highest soil δ 15 N, where α diversity was lowered. In a large scale seabirds may increase β or γ plants diversity by favouring the occurrence and enhancing the abundance of some nitrophilous plant species that are relatively widespread, although locally in areas that are not ornithogenically-influenced can have very low abundance on like. O. digyna, C. groenladica or C. arcticum.
The present study is of particular importance given the recent strong climatic and oceanographic changes in the Arctic which influence both marine and terrestrial parts of the ecosystem. Forecasts for zooplankton communities suggests these changes will lead to decreasing availability of Calanus glacialis-the staple food of little auks, and consequently a decline of the auk population [21,68,69]. This may lead to reduction in ornithogenic nutrient transfer onto land, impacting the tundra plant and animal communities [70,71] Supporting Information