Increasing atmospheric CO 2 concentrations outweighs effects of stand density in determining growth and water use ef ﬁ ciency in Pinus ponderosa of the semi-arid grasslands of Nebraska (U.S.A.) Global Ecology and Conservation

This study investigated the impacts of environmental (e.g., climate and CO 2 level) and ecological (e.g., stand density) factors on the long-term growth and physiology of ponderosa pine ( Pinus ponderosa ) in a semi-arid north American grassland. We hypothesized that ponderosa pine long-term growth patterns were positively in ﬂ uenced by an increase in atmospheric CO 2 concentrations and a decrease in stand density. To test this hypothesis, comparison of long-term trends in tree-ring width and carbon and oxygen stable isotopic composition of trees growing in dense and sparse forest stands were carried out at two sites located in the Nebraska National Forest. Results indicated that tree-ring growth increased over time, more at the sparse than at the dense stands. In addition, the carbon and oxygen isotopic ratios showed long-term increases in intrinsic water use ef ﬁ ciency (WUE i ), with little difference between dense and sparse stands. We found a clear trend over time in ponderosa pine tree growth and WUE i , mechanistically linked to long-term changes in global CO 2 concentration. The study also highlighted that global factors tend to outweigh local effects of stand density in determining long-term trends in ponderosa pine growth. Finally, we discuss the implications of these results for woody encroachment into grasslands of Nebraska and we underlined how the use of long-term time series is crucial for understanding those ecosystems and to guarantee their conservation. of including Pinus ponderosa . in site consists of trees that established over time in adjacent open grasslands. selected dense site thinning treatments, some tree mortality and self-thinning have and


Introduction
Spatiotemporal changes in grassland ecosystems processes are strongly influenced by natural disturbances, climatic variability and human activities, such as grazing management. Disentangling the role played by these factors on grassland ecological processes represents one of the major challenges of scientific community and forest managers. Changes in ecosystem processes are occurring at rates that are unprecedented in human history (Collins et al., 2011) and cover many broad geographic regions (Breshears et al., 2005;Goetz et al., 2012;Vayreda et al., 2012). Uncertainties on the trajectory of these changes under changing climatic conditions are particularly critical in a semi-arid environment where temperature, precipitation, and disturbances have long been viewed as key factors determining vegetation type, including transitions between grasslands, savannas and forests (Bond et al., 2005).
In Nebraska, as well as in other semi-arid regions of the Great Plains, one of the consequences of these changes is a shift in vegetation cover, such as the transformation of savanna-like ponderosa pine (Pinus ponderosa) ecosystems into dense forests, and its expansion into adjacent grasslands from historical grasslands-woodlands ecotones (e.g., Eggemeyer et al., 2006Eggemeyer et al., , 2009Msanne et al., 2017). Ponderosa pine (Pinus ponderosa P. & C. Lawson) is a tree species with one of the largest distribution area in North America (Sala et al., 2005), and recently it has expanded into grasslands from its original forest sites. Factors driving this phenomenon are still debated. Although, several factors have been indicated as drivers of those changes, such as decreased fire frequency (Bond et al., 2005), over-grazing, climate change, land use change, and human and natural enhanced seed dispersal (Awada et al., 2013 and references therein), it is difficult to disentangle such complex interactions between biotic and abiotic factors. Some studies have suggested that increasing atmospheric carbon dioxide (CO 2 ) concentrations may influence photosynthesis rates (Ausennac, 2002) and foster woody plant encroachment into grasslands (Stevens et al., 2016). Indeed, elevated CO 2 concentrations can enhance intrinsic water-use efficiency (WUE i ) of trees by reducing stomatal conductance (Polley et al., 1997;Bond and Midgley, 2000;Saurer et al., 2004;Battipaglia et al., 2013;Guerrieri et al., 2020) or increasing photosynthetic rates, which allow plants to increase their drought resistance (Drake et al., 2017). However, the role of increasing CO 2 concentrations, especially increasing plant WUE i and growth, in the physiological adaptation of plants to newly established climatic conditions in semi-arid ecosystems remains poorly understood (Archer et al., 1995;Davis et al., 1999;Dickie et al., 2007;Polley et al., 2003;Van Auken, 2009;Classen et al., 2010;Sullivan et al., 2017).
Tree-ring stable isotope ratios give insight into the ecophysiological processes triggering the trees response to past environmental conditions, including increase in CO 2 concentration in the atmosphere (Scheidegger et al., 2000;Barbour, 2007;Battipaglia et al., 2010;Sullivan et al., 2017). Once allocation patterns of carbon in trees are identified , carbon stable isotope ratio (d 13 C) can then be used to estimate annual variability of WUE i , since d 13 C values in organic matter reflect the ratio between the partial pressure of CO 2 in leaf intracellular space (c i ) and the partial pressure of CO 2 in the atmosphere (c a ) (Farquhar et al., 1982). Factors like limited water availability and increase in CO 2 concentrations in the atmosphere reduce stomatal conductance, triggering changes in WUE i and in c i that are recorded by variations in d 13 C in assimilated CO 2 and plant tissues (Seibt et al., 2008;Maseyk et al., 2011). The d 18 O of plant tissue depends primarily on source water signature and the ratio of ambient to leaf internal spaces, on evaporative enrichment at the leaf level and the following post-photosynthetic exchange rates during synthesis of carbohydrates (Roden and Ehleringer, 2000;Treydte et al., 2014). Changes in stomatal conductance (g s ) affect the evaporative enrichment of leaf water (Barbour, 2007) leaving a strong signal in the oxygen isotope ratio (d 18 O) values of tree rings. The combination of d 13 C and d 18 O provides complementary information, as d 18 O is linked to g s but it is not affected by net photosynthesis (A), and can thus allow the separation of the independent effects of A and g s on d 13 C and WUE i (Scheidegger et al., 2000;Roden and Farquhar, 2012;Battipaglia et al., 2013). Even if the interpretation of the double model d 13 C -d 18 O is not straightforward and may be sometimes hampered by source water isotope changes , it can still provide important information when applied in water-limited ecosystems (Moreno-Guti errez et al., 2012;Gessler et al., 2014;Altieri et al., 2015;Battipaglia et al., 2016).
In this study, we investigated the impact of stand density, climate variability, and increasing atmospheric CO 2 concentration on long-term growth and physiological processes of ponderosa pine trees in the semi-arid Nebraska grasslands, using tree-ring stable isotopes. We selected two types of stands with sharply contrasting structure and density: a dense ponderosa pine plantation, and a neighboring open woodland with scattered pines. We hypothesized that ponderosa pine growth at both sites is influenced by stand density and that an increase in atmospheric CO 2 concentrations would enhance tree water use efficiency, contributing to the ability of this species to cope with drought stress. Further, we discuss how elevated CO 2 could influence woody encroachment in such environment.

Study area
The study sites were located in the Nebraska National Forest (NNF), Halsey, Nebraska (825 m elevation, lat. 41 51 0 45 00 N, long. 100 22 0 06 00 W). The NNF is a 25,000 ha experimental forest, established in 1938 in the semi-arid C 4 dominated grasslands of the Nebraska Sandhills (McEntee, 1941). Tree planting started in 1920s and continued in subsequent decades. A total of~10,000 ha was hand planted with monocultures of coniferous species including Pinus ponderosa. Trees in sparse site consists mostly of volunteer trees that established over time in adjacent open grasslands. The selected dense site has not undergone thinning treatments, but some tree mortality and self-thinning have occurred over time. Temperature and precipitation data for the study period (1938e2014) for the two sites were obtained from the HPRCC (High Plains Regional Climate Center, HPRCC, University of Nebraska e Lincoln, http://www.hprcc.unl.edu, from local weather stations). The stations Halsey 2 W and Bessey Nebraska are in NNF, between 2 and 10 km North of the NNF sites, the station Dunning 6 NW 15 kme18 km Northeast of NNF sites. NNF climate stations Halsey 2 W and Bessey Nebraska are within 0.6 km of each other, whereas the Dunning 6 NW station is 13.5 km to the East (http://hprcc.unl.edu/).
The climate is semiarid, the mean annual temperature, recorded in the period 1950e2015, is 8.4 C, mean minimum temperature in January is À13 C and mean maximum temperature in July is 32 C. Mean annual precipitation is 570 mm. Vegetation season occurs from April to September, when around 75 percent of precipitation falls as rain. Soils are sandy (mixed, mesic, Typic Ustipsamments) (Lewis and Kuzila, 1998;Sherfey et al., 1965). The CO 2 records used for this study were derived from the ice cores obtained at Law Dome, East Antarctica, see Battipaglia et al. (2015) for details. The data are global but it has been demonstrated that CO 2 mixes well throughout the atmosphere, consequently, local CO 2 trend is statistically indistinguishable from the trend in global CO 2 levels (Thoning et al., 2015).

Sampling and tree-ring analyses
Three cores, from each of 20 ponderosa pine trees selected at each site, were taken with a 0.5 cm increment borer (Hagl€ ofs, Sweden) and with an angle of 90 between them. Two cores were used for ring-width measurements and a third one for isotopic analysis.
Tree-ring width (TRW) was measured at a resolution of 0.01 mm, using LINTAB equipment coupled with a stereoscope and with TSAP software (Frank Rinn, Heidelberg, Germany). After visual cross-dating (Schweingruber, 1996), samples from each site were cross-dated using the Gleichl€ aufigkeit, a statistical measure of the year-to-year agreement between the interval trends of the chronologies (Kaennel and Schweingruber, 1995). Student t-test was used to determine if the means of two curves are significantly different from each other.
The quality of the cross-dating and potential errors were estimated with the program COFECHA (Holmes, 1993). Corrected series were detrended with a cubic smoothing spline with a 50% frequency response cut-off of 10 years to remove long-term growth trends embedded in the raw ring-width series, which are induced by non-climatic influences, such as aging and competition between trees (Fritts, 1976). Ring-width indices were calculated as residuals from the estimated age trend. For each site chronology, the following parameters were calculated, using 30-year common interval with an overlap of 15 years: the signal-to-noise ratio (SNR), indicating the proportion of explainable variation divided by the unexplainable variation; the expressed population signal (EPS), indicating how the constructed chronology represented the hypothetical perfect population chronology, and the mean inter-series correlation (RBAR), measuring the common variance between the single series in a chronology (Table 1). Further, the mean sensitivity (MS), indicating the co-variation of the time-series and the series intercorrelation (S.I.), indicating the stand level signal, were reported (Table 1). We calculated the annual basal area increment (BAI) for each chronology, since it is less dependent on age (Biondi and Qeadan, 2008) and can be used to calculate cumulative basal area increment of each stand along the whole chronologies (Klein et al., 2014;Battipaglia et al., 2015).

Carbon and oxygen isotopes and water use efficiency
The 20 years prior to the date of sampling were investigated with annual resolution. Each annual ring was split and ground with a centrifugal mill (ZM 1000, Retsch, Germany) using a mesh size of 0.5 mm to assure homogeneity. Cellulose was extracted following the procedure described in Boettger et al. (2007) and Battipaglia et al. (2008). The carbon (C) and oxygen (O) stable isotopic composition was measured by continuous-flow isotope ratio mass spectrometry (Delta V Advantage, Thermo Scientific), the carbon by combustion in a EURO EA elemental analyzer (EuroVector, Milan, Italy), the oxygen by pyrolysis in a high-temperature oxygen analyzer (HEKAtech, Wegberg, Germany). The d 13 C raw data were corrected for the changing d 13 C of atmospheric CO 2 since the beginning of industrialization (Francey et al., 1999;McCarroll and Loader, 2004).
The corrected series were used in statistical analyses. The values of d 13 C in tree rings were used to calculate the intrinsic water use efficiency (WUE i ), defined as the ratio between net photosynthesis (A) and stomatal conductance (g s ) following the approach presented in Battipaglia et al., (2016).

Relationship between environmental conditions and tree growth
Differences of tree-ring width as well as of isotope composition in trees growing at sparse and dense sites were tested using a Student's t-test. All data were checked for normality before applying inferential statistics. For the whole period, in order to understand the effects of climate and increasing CO 2 , onto tree growth and physiology, we used a multivariate statistic known as Two-Block Partial Least Squares (2 B-PLS). This statistical approach can be successfully applied to matrixes with a low sample size and with highly correlated variables (Barker and Rayens, 2003;Carrascal et al., 2009), and, although originally designed for morphometric analyses, it was recently applied in several ecological contexts (Innangi et al., 2017(Innangi et al., , 2018(Innangi et al., , 2019. In contrast to other well-established multivariate statistics (e.g. Canonical Correlation Analysis), 2 B-PLS aims at finding latent variables that can explain the covariance between two multivariate matrixes, returning variables that account as much as possible for the covariation between two sets of variables (Rohlf and Corti, 2000).
The patterns of covariance between the two matrixes can be represented by a scatterplot for the first axis of the 2 B-PLS where the x-axis and the y-axis represent the two multivariate matrixes, respectively. These patterns of covariance can be interpreted by examining the correlations of the two multivariate matrixes (blocks) with first axis of the 2 B-PLS. Thus, patterns of positive or inverse correlation can be asserted both within and between matrixes, as variables that are positively correlated within one matrix and have the same sign on the other matrix highlight a pattern of covariance. Thus, the higher the correlation coefficient (regardless of its sign), the stronger the variable drives the scatterplot for the multivariate dataset within and between the considered blocks. All analyses were done in R 3.6.1 (R Core Team, 2019), using packages 'plsdepot' (Sanchez and Sanchez, 2012) and 'ggplot2' (Wickham, 2016).

Tree growth
Trees at the sparse and dense sites belonged to the same age class, with a mean of 80 ± 16 years. Site chronologies presented a high EPS values (>0.85) indicating that they can be considered representative of the radial growth of the studied stands (Wigley et al., 1984). Statistics (Table 1) revealed a strong common growth signal between the two sites and similar trends in terms of inter-annual changes linked to climatic fluctuations.
The tree ring-width (TRW) chronologies of sparse and dense sites presented similar trends despite the different stand densities. Further, significant correlation was found between stands (r ¼ 0.67, P < 0.01) during the whole growing period (Fig. 1) and in particular during 2000e2014 years (r ¼ 0.83, P < 0.001) when an increase of growth was recorded. Cumulative basal area increment (Fig. 2) for ponderosa pine showed a significant difference between dense and sparse site, and at both sites a statistically and highly significant increase over the years, especially at the sparse site (P < 0.001).

WUEi and d 18 O
The water use efficiency (WUE i ) of ponderosa pine, derived from tree-ring d 13 C measurements (Fig. 3), showed similar mean ± standard deviation values at the two sites (WUE i average sparse ¼ 112 ± 7 mmol mol À1 ; WUE i average dense ¼ 113 ± 7 mmol mol À1 ). A clear increase (P < 0.01) in WUE i was observed around the period 2000e2006 (WUE i average sparse ¼ 117 ± 5 mmol mol À1 ; WUE i average dense ¼ 119 ± 4 mmol mol À1 ), followed by a non-significant decrease. WUE i of sparse site had a positive correlation with temperature (r 2 ¼ 0.45, p < 0.05), with high WUE i during years characterized by high temperature (Fig. 3) At both sites a significant positive correlation between d 18 O and WUE i was found (r ¼ 0.25 in sparse site; p < 0,05 and r ¼ 020, p < 0.05 in dense site). The results from the first axis of the 2 B-PLS analysis are reported in Fig. 5. The x-axis of the plot represents the dendroecological data (Block 1, tree ring width, d 13 C, d 18 O) while Block 2, i.e., y-axis of the plot, represents environmental data (i.e. total precipitation, mean annual temperature and CO 2 concentration). A clear gradient from the lower left to the upper right quadrant of the plot can be observed. In detail, an array of data (between 1994 and 1999) can be seen in the lower left quadrant. This data demonstrates both higher TRW (in both sparse and dense stands) and an increase in total precipitation, as both TRW and P have negative signs for their correlation with axis 1 of the 2 B-PLS, thus being positively correlated between datasets. Data from years 2000e2014, are grouped in the upper right quadrant. In general, these data were on the opposite quadrant from 1994 to 1999 data, implying a decrease in both TRW and total precipitation. In addition, these data were characterized by a simultaneous increase in CO 2 and d 13 C and d 18 O, although with some important distinctions. The covariance between CO 2 and d 13 C was sharper in the dense stands, while d 18 O values in the sparse stands were more prominent. Mean annual temperature and d 13 C and d 18 O at the sparse and dense stands showed little or no impact on growth. 4. Discussion

Effect of environmental conditions on ponderosa pine growth
In this study, we found an increase in ponderosa pine growth starting from 1960 at both sites, and in particular at the sparse one. These results agree with previous studies suggesting that tree species' responses to climate are site-dependent and may be modulated by stand density (Moreno-Gutierrez et al., 2012;Camarero et al., 2013;Sanchez-Salguero et al., 2013). Zalloni et al. (2019) found that intraspecific facilitation effects decreased with increasing stand density in a Quercus ilex stands, while Niccoli et al. (2020) demonstrated that sparse stands, resulting from selective thinning, have higher productivity associated with increased photosynthetic rate and decreased water loss relative to low thinned stands.  Moreover, several microclimatic parameters such as soil nutrient availability were shown to change with increased stand density (Barron-Gafford et al., 2003;Zhao et al., 2012;Liu et al., 2013;Bolat, 2014), impacting tree growth at site level (Granda et al., 2017). In areas with poor availability of nutrients in the soil, or low N concentrations, as in the Nebraska sites (Awada et al., 2013), a low stand density may be enhancing tree-growth rates (Fern andez-Martínez and Fleck, 2016).
When considering the environmental factors triggering tree growth, neither precipitation nor temperature seemed to be strongly influencing, while the most important driver of tree growth appeared to be CO 2 concentrations. There are still contrasting evidences on the role of elevated CO 2 concentrations on tree growth. Some studies report an increase in wood biomass in forest ecosystems worldwide (Yu et al., 2019), whereas many others report no fertilization effects on tree growth (Peneuelas et al., 2011;Battipaglia et al., 2015, van der Sleen et al., 2015Brienen et al., 2016).
However, recently, it has been hypothesized that a global factor, such as atmospheric CO 2 enrichment could be considered an important triggering factor for tree growth overriding the local environment, especially when trees coexist with C 4 grassland species, as in our study sites, the so-called encroachment process (Stevens et al., 2016 and reference therein; Sullivan et al., 2017). Indeed, it has been proposed that elevated atmospheric concentration of CO 2 favors C 3 in comparison to C 4 plants through two main processes: i) increasing of carbon allocation to roots or ii) increasing soil water availability due to reduced plant transpiration (Devine et al., 2017), potentially favoring trees over grasses. Furthermore, young trees in open environments, such as at the sparse sites, grow in relatively resource-rich environments relatively to trees growing in denser sites where competition for water and light is high and may therefore be expected to benefit the most from CO 2 fertilization (Bond and Midgley, 2000;Msanne et al., 2017), especially in water-limited sites as studied sites (Aus der Au et al., 2018).

Carbon uptake and water use efficiency
Ponderosa pine trees growing in Nebraska National Forest showed an increase in WUE i that was significant between 2000 and 2006, during a period of several years of consecutive drought, (High Plains Regional Climate Center, University of Nebraska-Lincoln, https://hprcc.unl.edu/) and among the environmental drivers of d 13 C and 13 C-derived WUE i change. The increase in atmospheric CO 2 concentrations proved to be the most relevant in our study site. High temperature influenced mostly the increase of WUE i in the sparse site and only in some years (such as 2002, 2006, 2012). Further, during the investigated period in the study, no change in land use or fire regimes was reported (https://www.fs.usda.gov), or variation of other environmental factors (Aus der Au et al., 2018). The dual-isotope approach (Scheidegger et al., 2000) has commonly been applied in tree-ring studies to trace the balance between carbon uptake and water losses (e.g., Moreno-Gutierrez et al., 2015). d 13 C of tree rings depends on factors affecting CO 2 assimilation, mainly controlled by photosynthetic rate (A) and stomatal conductance (gs), expressed by A/gs (intrinsic water-use efficiency; Farquhar et al., 1989). On the other hand, d 18 O can be considered mainly an indicator of stable isotope composition of source water, and leaf water enrichment due to stomata transpiration (Labuhn et al., 2014;Treydte et al., 2014). In fact, both isotopes are influenced by changes in stomata conductance, but because of the fact that d 18 O is not directly linked to photosynthetic activity, the combination of d 13 C and d 18 O helps clarify the role of stomatal activity and photosynthesis on tree performance, which may ultimately determine changes in productivity (Barbour et al., 2002;Barnard et al., 2012;Shestakova et al., 2017). The d 13 C and consequentially the 13 C-derived WUE i increase, reported in our study during the last years, may be caused by an enrichment in A, given that d 18 O-derived gs was maintained rather stable for the whole study period at both sites (Scheidegger et al., 2000), as also shown by the strong observed basal area increment-increase. Similarly, Fern andez-Martínez and Fleck (2016) found increased photosynthetic rates and WUE i in Pinus uncinata Mill without changes in gs under elevated CO 2 . This could indicate that ponderosa pine trees can benefit of higher CO 2 concentrations in the atmosphere and warm temperatures once established, enhancing A without suffering water shortage. It has been demonstrated that ponderosa pine possesses avoidance strategies that allow it to physiologically tolerate soil moisture stress in semiarid grassland environments (Martinez-Vilalta et al., 2004;Eggemeyer et al., 2006). The deep roots of ponderosa pine allow this species to access water to 12 m deep in the soil profile (Burns and Honkala, 1990) below the frozen layer during winter (Eggemeyer et al., 2009), and to avoid drought stress during the summer (Eggemeyer et al., 2009), thus facilitating its survival in semi-arid environments (Kolb and Robberecht, 1996). The positive correlations found between WUE i and growth, unlike those found in other studies (e.g. Penuelas et al., 2011), indicate that the increase in atmospheric CO 2 concentrations led to increased C inputs in both sites. We found that the two sites were different in terms of growth with the sparse possessing higher growth rate relative to the dense site. However, no significant differences in WUE i or d 18 O were observed between the two sites, and the dense stand did not show any signs of major stress linked to the intense inter-tree competition for soil moisture or significantly higher canopy interception of precipitation (Chirino et al., 2006). Previous study by Eggemeyer et al. (2009) reported the importance of root plasticity in this species, with hydrogen and oxygen isotopes showing that the species competes with grasses for water in the top soil (0e30 cm) during spring (period of high precipitation, May and June) and shifts its water uptake to deeper layers of the soil profile during July and August (period of drought). On the contrary in Mediterranean environments, Linares et al. (2010) found that trees of Abies pinsapo were suffering from intense inter-tree competition (a long-term stress), and Moreno-Gutierrez et al. (2012) found in two Aleppo pines (Pinus halepensis) stands that tree-ring growth and d 18 O, but not d 13 C, were affected by stand structure in severely water-limited ecosystem. Similar d 13 C values between dense and open sparse stands suggest that ponderosa pine maintains in the long term a homeostatic control of the ratio ci/ca (the intercellular to atmospheric CO 2 concentration, which determines d 13 C and WUE i ) as a constant uptake of water (Eggemeyer et al., 2006;McDowell et al., 2006;Maseyk et al., 2011). These results are based on the "linear" carbon isotope plant fractionation model, but more accurate results might be achieved with the detailed fractionation model including mesophyll conductance that is also influenced by environment (e.g., Seibt et al., 2008). However, those models are difficult to estimate in hindsight and if unaccounted for, they could induce some additional uncertainty. A further bias could be linked to carbon and oxygen signals recorded in tree rings. Indeed, more than one carbon source is integrated in tree-ring d 13 C, from leaf photoassimilates to stored C pools in wood (Seibt et al., 2008). Further, the use of different source waters along the growing season (Eggemeyer et al., 2009) could influence the signature of d 18 O in tree rings (Roden and Siegwolf, 2000).

Linking growth and physiological processes to pine encroachment
As woody species invade and alter species composition, the source and amount of water uptake by plants can be altered, affecting not only the site water balance but also the water available to native trees and grasses (Awada et al., 2013(Awada et al., , 2019. This allows the woody species to better survive on semi-arid sites and to uptake water sources that usually are beyond the reach of grass vegetation. A recent study on 115-year long selection of grasses and woody species collected in new Mexico, showed that elevated CO 2 influenced plant sensitivity to water shortage, increasing tree WUE i in arid environment and making them less vulnerable to water stress (Drake et al., 2017).
Our findings indicated the crucial role of elevated CO 2 concentrations in maintaining ponderosa pine at the semi-arid edges of its range and helping it to save water without paying huge carbon penalty during photosynthesis, increasing its growth and productivity. Those findings open important and timely questions on the future conservation of those ecosystems where encroachment processes could be triggered by climate change and in particular by increasing CO 2 concentrations (Van Auken, 2000; Morgan et al., 2007;Staver et al., 2011;O'Connor et al., 2014).

Authors contributions
PC, TA, GB conceived and designed the study. RADA and GB performed sampling and analyses. PC, TA, MS, MI contributed to the analysis tools. All authors contributed to interpretation of the overall data. GB wrote the main part of the manuscript. All authors approved the submitted version of the manuscript.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.