Nitrogen cycling in an extreme hyperarid environment inferred from δ15N analyses of plants, soils and herbivore diet

Climate controls on the nitrogen cycle are suggested by the negative correlation between precipitation and δ15N values across different ecosystems. For arid ecosystems this is unclear, as water limitation among other factors can confound this relationship. We measured herbivore feces, foliar and soil δ15N and δ13C values and chemically characterized soils (pH and elemental composition) along an elevational/climatic gradient in the Atacama Desert, northern Chile. Although very positive δ15N values span the entire gradient, soil δ15N values show a positive correlation with aridity as expected. In contrast, foliar δ15N values and herbivore feces show a hump-shaped relationship with elevation, suggesting that plants are using a different N source, possibly of biotic origin. Thus at the extreme limits of plant life, biotic interactions may be just as important as abiotic processes, such as climate in explaining ecosystem δ15N values.

. Study site, climate and sampling sites in northern Chile. Regional context of northern Chile showing location of the Salar de Atacama and adjacent Andes (right inset) and a digital elevation model indicating where our sampling sites (lower inset, colored dots) are along the Talabre-Lejía Transect (TLT). Isohyets are in the same corresponding color as the respective sampling sites. The software used to create the map was QGIS 2.10 with Openlayer plugin, STRM30 53,54 elevation model (Data: SIO, NOAA, U.S. Navy, NGA, GEBCO) and Landsat 8 Satellite image (Data available from the U.S. Geological Survey).
We surveyed and collected plants, surface soil samples and herbivore feces to assess impacts on diet along the 1800 m elevational gradient (TLT) that ranges from the extreme hyperarid conditions at 2700 meters above sea level (m) to extremely cold conditions at 4500 m. We obtained soil physical (%sand, %clay, %silt) and chemical properties, including total N, NO 3 , NH 4 , pH, organic matter (OM), macronutrients (P, C, S, K), SAR (sodium adsorption ratio) and analyzed our soil samples for δ 15 N to better understand the factors that could be behind N cycling as well as soil development along this climatic/altitudinal gradient.
We describe the relationship between soils, foliar samples and herbivore feces with climate across this gradient and construct the isotope landscape for this extreme environment. We evaluate the mechanisms proposed to explain these relationships, divided into those that focus on abiotic factors (i.e. climate, pH) and those that focus on biotic factors (i.e. diversity, presence of N 2 -fixers and how these relate to microorganisms). Finally, we compare soil, plant and herbivore feces δ 15 N values to infer how N acquisition occurs in plants from these extreme environments.

Results
Plant and soil descriptions. MAP, MAT and therefore aridity are mainly a function of elevation in this area of the Atacama. All of these variables co-vary and we were unable to discern which had the most explanatory power. Hence, we present our results in relation to elevation in the understanding that this co-varies positively with MAP and aridity (R 2 = 0.93, p < 0.001) and negatively co-varies with MAT (R 2 = − 0.99, p < 0.001). Plant species richness and percent relative cover exhibit a "hump-shaped" curve with elevation (Fig. 2). The calculated De Martonne aridity index (0-5 hyperarid, 5-10 arid and 10-20 semiarid) 30 varied between 11.4 at 4500 m and 0.5 at 2700 m (Table 1). Zones characterized as hyperarid by this index are approximately those between 2700 and 3900 m.
A pronounced pH gradient occurs along the TLT, with acidic soils occurring at high elevations (pH 5.23) and alkaline soils at lower elevations (pH 8.54) ( Table 1). Soil pH shows significant correlations with δ 15 N, δ 13 C, aridity, mean C (mg/kg) and some macronutrients (N, P, S, K) ( Table 2). The soil physical composition is mainly sand ~81%, followed by silt ~12% and finally, clay ~7%. A negative correlation (R 2 = − 0.66, p < 0.005) occurs between % clay and the aridity index (a higher index means less aridity). Total nitrogen was low throughout the environmental gradient (0.26-0.60 mg/g) and shows no significant correlation with aridity (R 2 = 0.19, p > 0.05) ( Table 2). Other soil parameters (total N, NO 3 , C and P) correlate with biotic factors and do not correlate with aridity. These show increases at mid-elevations, where temperature is not as cold as at higher elevations and rainfall is not as dry as at lower elevations. Hence, plant diversity, plant cover and organic matter (OM) are highest at mid-elevations along our gradient (Table 2 and Fig. 2). P increases with MAP, whereas K decreases. A high sodium adsorption ratio (SAR) implies poorly irrigated soils and as expected, correlates positively with aridity. Soil isotopes. The isotopic δ 15 N soil values range from 3.3% to 12.2% with a δ 15 N soil mean of 8.2% (n = 40) ( Table 3 shows the mean site δ 15 N soil values, the entire list of 40 surficial soil samples can be found in Supplementary Table 1). Regression analysis shows that δ 15 N soil values are strongly dependent on elevation, which explains 72% of the variance (Fig. 3, R 2 = 0.72, p < 0.001).
No apparent differences were observed between soil samples taken in different years after the rainy season (April 2011, 2012 and 2013) (Supplementary Table 1). We also compared how δ 15 N soil values changes with soil depth (0, 25 and 50 cm) at three representative sites (2900, 3600 and 4300 m) and found no apparent differences (Supplementary Table 3).
Herbivore feces isotopes. Herbivores feces, mainly rodents (N = 10) and camelids (N = 9), show a humpback relationship with elevation (R 2 = 0.67, p < 0.001) with no differential elevational responses across the different species sampled. The δ 15 N values (Table 5) ranges from 2.1% to 8.3% (Fig. 4). Plant δ 15 N values predict 53% (R 2 = 0.53, p < 0.005) of the variance observed in herbivore δ 15 N. Although plants and herbivores show similar signals (Fig. 4), at mid-elevations, herbivores (as expected by trophic enrichment) tend to be enriched in 15 N by ~1% compared to plants. At the extremes of our gradient, however, plants are enriched compared to herbivores, suggesting that herbivores at these elevations are mixing or selecting their plant sources from other local sources (such as a spring or deep canyon).

Discussion
Five different mechanisms have been proposed to explain how soil δ 15 N values relate to aridity or climate. These include: (1) high MAP increases soil organic matter which becomes depleted in 15 N in response to fractionation during mineralization 31 ; (2) increase in the discrimination of 15 N with increasing N soil reserves 17 ; (3) elevated temperatures in arid and hyper-arid environments can trigger increased N volatilization and the preferential removal of 14 N 32 ; (4) changes in the relative importance of within-system cycling versus inputs/outputs or the "openness" of the N cycle 1,13 ; and (5) the relationship between soil δ 15 N and climate is indirect and mediated through climatic effects on soil C concentration and clay amount 6 .
All of our sites exhibited positive δ 15 N soil values (1.9% to 11.6%%), indicating the preferential removal of 14 N from soils, coupled with overall lack of organic matter. The large ~10% variation is unusual and appears to be related to the extreme hyperaridity of our study system. As expected, the most positive values appear at the driest sites (~2700 m) and soil % N does not change across the TLT (Table 3). Because the organic matter and N content in our sites do not correlate with δ 15 N soil values (Table 2), the aforementioned mechanisms (1) and (2) are unlikely to explain our observations.
Mechanisms (3), (4) and (5) are more coherent with our results. Mechanism (3) could explain the very positive values found in the driest sites which may result from high N volatilization during elevated diurnal temperatures and could be associated with very alkaline soils (~pH = 9) 33 . However, TLT δ 15 N soil values may also increase due   by mechanism (4). Plant N demands increase with rainfall, creating a larger soil organic N pool. The retention of organic N, instead of inorganic labile pools, reduces the N outputs and increases within-system cycling (i.e. decreasing cycle "openness"). Although leaching would increase compared to the drier sites, organic N is less labile. Increased aridity, salinity and extreme pH decrease the proportional flux of ecosystem N into organic pools and increase the inorganic N pool and outputs. Since nearly all N transformation processes fractionate preferentially for 14 N 16 , cycle "openness" drives ecosystem N towards 15 N-enrichment (relative to atmospheric N 2 ) in the driest sites. Finally, % clay and carbon (mechanism 5) could explain in part some of the variability. Although soil C concentration exhibits little variation across the gradient (0.89 to 3.33 mg/kg), a significant correlation (R 2 = − 0.53, p < 0.05) with δ 15 N soil values does occur as expected. Percent clay in our soils varies from 3.3 to 11.9% and the correlation with δ 15 N soil values is significant (R 2 = 0.67, p < 0.005) yet the relation between MAP and % clay is actually inverse (R 2 = − 0.65, p < 0.005) as the amount of clay decreases with elevation. This implies that although % clay is a contributing factor, precipitation and the indirect effect as observed in global climate gradients 6 may not be the actual mechanism in our gradient from the Atacama. Furthermore, our soil δ 15 N samples exhibit stronger direct correlations with climatic variables (MAP and MAT) than with C or % clay (Table 2).   Other drivers besides climate have been proposed to explain δ 15 N soil variations including soil depth 34 , age 2 and the presence of N 2 fixers 35 , among others. Soil depth would not affect our results as all samples were taken from the first 5 cm and we further tested for variations in δ 15 N soil with depth (0, 25 and 50 cm) at three sites and no differences were found (Supplementary Table 3). A few N 2 fixers (Adesmia spp., Lupinus spp. and Hoffmannseggia doelli) are found along the TLT. These N 2 fixers are characterized by comparatively low δ 15 N values (close to atmospheric N 2 ~ 0%) 36 . Some of our sites were dominated by Lupinus spp. (L. subinflatus and L. oreophilus) during wet years, but the presence of these plants was not reflected in our δ 15 N soil values, which are highly enriched compared to atmospheric N 2 .
Ewing and collaborators 37 proposed a transformation in the N cycle from arid to hyperarid sites where soil N loss decreases to be almost negligible. They described two modes of N accumulation. In semiarid environments, total soil N increases with rainfall and is mostly organic. In arid to hyperarid sites total soil N increases with aridity and is mostly inorganic. We thus expected to find more N accumulated in the soils at our most arid sites. We also expected N to increase at the lower sites as the C/N ratio decreased and biological uptake of N ceased. Surprisingly, we found that although plant cover and total organic matter content decreases at the extremes of our survey, total soil N did not change despite a relative increase in total soil C with elevation (Table 3).   The mean soil δ 15 N value from all sites (8.2 ± 1.8%) concurs with what was expected for such an hyperarid area as extrapolated from a global dataset 6 . Within the gradient, we found an overall ~3% shift that correlates with an order of magnitude shift in MAP, which is more pronounced than the expected ~2% shift inferred from the global dataset 6 . We observed a linear relationship between elevation and soil δ 15 N values (R 2 = 0.72, p < 0.001), and a hump-shaped curve between elevation and foliar δ 15 N values (R 2 = 0.58, p < 0.001), (Figs 3 and 4). A hump-shaped second order polynomial is indeed a better fit than a first-order model to our foliar samples (p < 0.05).
These results imply that both biotic (e.g. plant cover, richness) and abiotic (climate gradients) factors constitute possible explanations for the isotopic properties observed. Indeed, δ 15 N foliar and δ 15 N herbivore exhibit a hump-shaped relationship with elevation (Fig. 4) which points to different drivers below and above ~3300 m. For sites > 3300 m, foliar δ 15 N appears to be coupled with soil δ 15 N, but < 3300 m (where vegetation cover falls to almost zero), these values are much less positive (~4%) than their corresponding soil δ 15 N values (Fig. 3). This apparent "decoupling" could be due to short-term plant δ 15 N dynamics of N cycling versus the long-term dynamics of soil δ 15 N values 38 .
To evaluate if the soil and plant values are effectively coupled or decoupled along the different parts of the gradient or across vegetation belts, we normalized and divided the data into three vegetation belts (prepuna, puna and steppe) (Supplementary Fig. 2). A t-test comparing soils and plants from each vegetation belt reveals no significant differences (p > 0.05) although plants fit a "hump-shaped" distribution better than our soil samples, which are clearly linear with elevation. Nevertheless, these results could be more robust as the number of data within each class is still low, an important prospect for future studies.
A recent study in arid grasslands from China 20 also described a hump-shaped curve between plant δ 15 N and aridity, similar to our results but with a different threshold. In that study 20 the tipping point occurs at sites with less than ~200 mm/yr of rainfall. Yet in our study, all sites are below 200 mm and the threshold is closer to 50 mm/yr (~3500 m). They point to a change in net plant N accumulation relative to gaseous N losses (volatilization) as the primary cause in determining the negative relationship between aridity and δ 15 N at the drier sites. The implication is that the net ecosystem N retention rate increases with aridity. This clearly was not the case for the Atacama as soil total N (mg/kg) does not correlate significantly with elevation, MAP or MAT (Table 2).
In animal tissues, δ 15 N values are derived from diet, climate and/or physiology, but there is some disagreement regarding which of these factors is (are) the main driver(s) in the observed 15 N abundance across environments or taxa [39][40][41] . It has also long been recognized that δ 15 N in animal tissues and feces are enriched compared to their diet (~3%) 39,42,43 . Our herbivore feces do not show any consistent 15 N enrichment compared to plants for those same elevations along our TLT (Fig. 4). In fact, feces 15 N abundance is overall very similar to the plants at the same elevations and displays a hump-shaped relationship when plotted against elevation, MAP, or MAT (Fig. 4). These results suggest that herbivores feces from the Atacama have an important component of non-digested plants, as seen by the presence of well-preserved fibers and plant remains in rodent feces 44 .
As expected, there is an overall strong negative correlation between the aridity index (which increases with decreasing "aridity") and soil δ 15 N values (R 2 = − 0.71, p < 0.001) ( Fig. 3 and Table 2) and an overall positive correlation in δ 15 N values between soil-plants-herbivores, as plants provide the only exogenous source of N to herbivores 41,45 . This relationship breaks down, however, between 3200 and 2700 m, where total vegetation cover drops to < 1%. Thus, plants and feces show an enrichment in 15 N with elevation from 2700 to 3500 m and a depletion in 15 N with elevation from 3500 to 4000 m (Fig. 4). These results suggest that different drivers for the N isotope natural abundance appear to be at work in the Atacama at different elevations and that feces δ 15 N values are directly established by diet and do not appear to be further enriched by animal metabolism. Previous research has also pointed out that in hyperarid environments, the typical positive relationship between aridity and δ 15 N may not be always accurate 41,46 .
We speculate that perhaps the role of microorganisms at such hyperarid sites in the Atacama could also be an important factor, either as free-living soil organisms or in association with plants 47 . Microorganisms are clearly important in fixing atmospheric N 2 and recycling organic N from leaf litter, processes that are known to deplete 15 N absorbed by plants and could explain the differences between soil δ 15 N and foliar δ 15 N at our lower elevation sites. Hence, plants could survive above micro-nutrient "islands" with different N sources compared to the surrounding soils in such an extreme environment.
Finally, the extreme positive values seen in our results can help shed further light on the mechanisms underlying the relationship between climate and δ 15 N. The soil δ 15 N pattern that occurs in our gradient appears to favor two different but complimentary hypotheses. As expected, soil δ 15 N values are correlated with elevation and climate, especially above 3300 m. As the driest desert of the world, water is the most important factor controlling plant cover and even diversity. The interplay between precipitation and temperature, however, regulates the presence of vegetation at higher elevations. In contrast, high diurnal temperatures and elevated pH at our most arid sites (~2700 m) would promote ammonium volatilization. Coupled with a more "openness" of the N cycle, this would further increase the loss of labile inorganic N pools through volatilization.
In contrast, the hump-shaped foliar and herbivore feces δ 15 N relationship with elevation points to different drivers along the gradient. Above 3300 m, foliar δ 15 N values are coupled with soil δ 15 N values and are likely regulated by climate and other abiotic factors. In contrast, below 3300 m foliar δ 15 N values become apparently "decoupled" from their corresponding soil δ 15 N values ( Supplementary Fig. 2) indicating biotic factors are more prevalent. Our results suggest that greater attention should be devoted to understanding the role of the microbiome in these extreme environments, where they could be important drivers of soil N availability.

Methods
We characterized 20 sites every 100 m of altitude ( Fig. 1 and Table 1) during four consecutive years (2011-2014), in April (after the rains, during the growing season). Maximum elevation was 4500 m (MAP ~ 160 mm) and the lowest elevation was at 2700 m (MAP ~ 10 mm) (Fig. 1). The De Martonne aridity index was calculated at each site as described in previous reports 30 (Table 1). Plant richness and plant cover were estimated using the McAuliffe log-series survey 48 with two plots of 250 m 2 for each of the 20 sites along the TLT (Fig. 2).
We analyzed two independently collected soil samples from each site. We collected each isotopic soil sample by removing the top 1 cm of surface soil and collecting the next five centimeters (from five different soils in each site). Samples were taken to the laboratory in sealed plastic bags, dried at 50 °C, and sieved and inspected under a binocular microscope to remove visible plant remains such as rootlets, before the δ 15 N and δ 13 C isotopic analyses. Table 3 shows mean values and the entire list of 40 soil samples can be found in Supplementary Table 1. Foliar samples were collected from various individual ramets of each species and homogenized for total C, N and δ 13 C, δ 15 N analyses (Table 4, mean values). As it is impossible to collect the same plant species across the entire gradient, six species were selected for isotope analyses (Baccharis tola, Jarava frigida, Maihueniopsis camachoi, Parastrephia quadrangularis, Tiquilia atacamensis and Atriplex imbricata) based on their widespread distribution within the Atacama gradient. The relative % plant cover of these species was at least 50% at each site. Leaves, fruits or seeds were sampled from 66 specimens and analyzed for total C, N and δ 13 C, δ 15 N (Supplementary Table 2). The plants were collected, pressed and dried in the field, and later placed in a drying oven at 50 °C in the laboratory (the entire list of 66 samples can be found as Supplementary Table 2). We also collected and analyzed 19 feces samples collected along the TLT during April 2011-2014 (Table 5). Rodents (Abrocoma and Phyllotis) are the most common taxa sampled (N = 10), but we also collected camelid feces (N = 9). δ 15 N was standardized with N 2 -Air and δ 13 C values to Vienna Pee Dee Belemnite (VPDB). The δ values were measured in units of per mil (%). A total of 36 soil and 50 foliar samples were submitted to the Cornell University Stable Isotope Laboratory (COIL) and analyzed on a Thermo Delta V Advantage isotope ratio mass spectrometer (IRMS) coupled with a NC2500 Elemental Analyzer. Herbivore feces (n = 19) and additional soil (14) and foliar (16) samples were measured at the Laboratory for Biogeochemistry and Applied Stable Isotopes (LABASI) of the Departamento de Ecología, Pontificia Universidad Católica de Chile using a Thermo Delta V Advantage IRMS coupled with a Flash2000 Elemental Analyzer (see Supplementary Tables 1 and 2). Overall, agreement between these two labs was tested by same sample comparisons and reproducibility was within the analytical error of these instruments (± 0.2%).
Soil samples for chemical analysis parameters were also collected from all stations along our survey during April 2012, 2013 and 2014. These data are used to calculate the correlation matrix (Table 2) as an average from these three years. Soil composition can change across the landscape, and also within sites, according to soil moisture, soil type, topography, aspect or vegetation, particularly in arid ecosystems, thus generating a heterogeneous biogeochemical landscape 49 . To handle this heterogeneity, the soil samples from each station consisted of the soil collected and mixed from 10 randomly placed quadrants (diameter: 15 cm, depth: 5 cm) (Supplementary Table 1).
Soil texture and chemical analysis was carried out by the Laboratorio Agroanalisis UC, Facultad de Agronomía e Ingeniería Forestal, Pontificia Universidad Católica de Chile, according to the methods established by the CNA of the Chilean Society of Soil Science 50 . Analyses included grain size, elemental composition, and soil pH, among others (see Table 2). Statistical Analysis. Simple linear models were used to study correlations and variation in δ 15 N and δ 13 C values. The soil and climatic parameters (pH, N (mg/g), C (mg/kg), MAP and MAT among others) were correlated with Pearson's correlation coefficient (Table 2). To analyze the coupling of the soil and foliar samples across the different vegetation belts or elevation we normalized each dataset, which are then graphed as box plots ( Supplementary Fig. 2). All analyses, charts and maps were performed in the R Programming Language and Quantum GIS software 51,52 .