Assessing the Contrasting Effects of the Exceptional 2015 Drought on the Carbon Dynamics in Two Norway Spruce Forest Ecosystems

: The occurrence of extreme drought poses a severe threat to forest ecosystems and reduces their capability to sequester carbon dioxide. This study analysed the impacts of a central European summer drought in 2015 on gross primary productivity (GPP) at two Norway spruce forest sites representing two contrasting climatic conditions—cold and humid climate at Bílý K˘rí˘z (CZ-BK1) vs. moderately warm and dry climate at Rájec (CZ-RAJ). The comparative analyses of GPP was based on a three-year eddy covariance dataset, where 2014 and 2016 represented years with normal conditions, while 2015 was characterized by dry conditions. A signiﬁcant decline in the forest GPP was found during the dry year of 2015, reaching 14% and 6% at CZ-BK1 and CZ-RAJ, respectively. The reduction in GPP coincided with high ecosystem respiration (R eco ) during the dry year period, especially during July and August, when several heat waves hit the region. Additional analyses of GPP decline during the dry year period suggested that a vapour pressure deﬁcit played a more important role than the soil volumetric water content at both investigated sites, highlighting the often neglected importance of considering the species hydraulic strategy (isohydric vs. anisohydric) in drought impact assessments. The study indicates the high vulnerability of the Norway spruce forest to drought stress, especially at sites with precipitation equal or smaller than the atmospheric evaporative demand. Since central Europe is currently experiencing large-scale dieback of Norway spruce forests in lowlands and uplands (such as for CZ-RAJ conditions), the ﬁndings of this study may help to quantitatively assess the fate of these widespread cultures under future climate projections, and may help to delimitate the areas of their sustainable production.


Introduction
Terrestrial ecosystems, such as forests, play a significant role in the global carbon cycle by sequestering carbon dioxide (CO 2 ) through photosynthesis and the conversion of biomass into stable soil organic compounds [1][2][3][4][5][6]. However, carbon uptake and carbon allocation can be reduced by high vapour pressure deficits (VPDs) and soil water deficits, commonly observed during periods with increased temperature and lack of precipitation [7,8]. These physiological responses to drought in plants vary at both local and global scales, depending on the type of plant species (based on the different levels of resilience to water stress), the local climatic condition, and other additional factors [9][10][11][12][13]. It has also been widely demonstrated that climate change increases the likelihood and severity of such drought events, especially within the temperate regions [14][15][16][17][18]. Thus, despite the resilience of some tree species, the projected increase in the occurrence and severity of extreme drought events will further reduce forest ecosystem photosynthesis in many parts of Europe, especially at unsuitable locations with sub-optimal climatic conditions [19].
Drought does not have a single definition and can be used in different contexts. However, typically, it is used to describe the periods with high VPD (meteorological drought) and low soil moisture (edaphic drought) [20]. The severity of the drought impact on carbon exchange depends on the site characteristics, duration and intensity of the drought periods [18,21,22]. In Europe, the 2003 summer drought period resulted in an estimated loss of about 30% in GPP (0.5 P gC y −1 ) over both the northern Mediterranean forests and temperate deciduous beech forests [23,24]. This was the result of the prolonged abnormal reduction in the soil moisture and high air temperatures below the wilting point and high air temperatures recorded across Europe during the summer period in 2003. The summer drought of 2015 was considered one of the most severe drought events in Europe after the 2003 summer drought [25]. During this period, much of the European continent (Poland, the Czech Republic, Slovakia, Western Ukraine, and Belarus) was severely affected in June and July 2015, due to the exceptionally high maximum daily air temperatures and the prolonged rainfall shortage by 31% since April [26][27][28][29]. Consequently, this had a significant effect on forest growth, as observed in both spruce and beech forests, within the Czech Republic. Nonetheless, spruce forest showed a higher sensitivity to drought as compared to beech forests [30,31].
The Norway spruce (Picea abies (L.) H. Karst) has significant economic and ecological importance within Europe [32]. However, its shallow root system makes it vulnerable to drought stress, especially at lower altitudes [33]. Since spruce thrives well in cold and humid regions, a significant decrease in soil moisture coupled with high air temperatures could have severe consequences on the spruce forest growth [34]. Therefore, the 2015 summer drought episode provided a good opportunity to further examine the response of spruce forest stands with different climatic conditions to extreme drought stress events.
The eddy covariance (EC) technique provides direct measurements of CO 2 exchange between the atmosphere and the underlying ecosystem [35]. It is a convenient and widely used approach for observation of a forest stand carbon uptake and its dynamic response to environmental variables. Obtained CO 2 fluxes are integrated over hundreds of square meters and resolved to half-hourly intervals [35,36]. Continuous micrometeorological and EC measurements over multiple years allow to detect drought periods and evaluate their impacts on CO 2 exchange.
In this study, we hypothesized that the impact of extreme drought conditions will be more severe on the Norway spruce forest stand located in dry and hot climates, due to higher VPD and the sub-optimal supply of soil moisture in such drier climates. To test this hypothesis, we sought to evaluate the different effects of the 2015 drought during the main growing season in the wet and dry spruce forest ecosystems, compared to normal climatic conditions within the same period of two adjacent years (2014 and 2016). Secondly, we aimed to determine the influence of critical site-specific environmental factors (such as the VPD and soil volumetric water content) on the drought stress response at each of the spruce forest stands located in mountainous (around 900 m a.s.l) and highland (up to 600 m a.s.l) regions within the Czech Republic, using long-term eddy covariance CO 2 flux data.

Station Description
The study used multiyear (2014-2016) measurements from Bílý Kȓíz and Rájec ecosystem stations that are part of the CzeCOS (Czech Carbon Observation System; http:// www.czecos.cz/ (accessed on 6 March 2021)) and FLUXNET (Flux Tower Network; https: //fluxnet.fluxdata.org/ (accessed on 6 March 2021)) station network. The FLUXNET station IDs for the wet and dry spruce forest are CZ-BK1 and CZ-RAJ, respectively. CZ-BK1 is also a candidate for ICOS (Integrated Carbon Observation System; https://www.icos-cp.eu/ (accessed on 6 March 2021)) network. The main vegetation cover in both stations is the even-aged stand of Norway spruce. Their main characteristics are presented in Table 1. CZ-BK1 is situated in the Moravian-Silesian Beskids Mountains in the Czech Republic and has a cold and humid climate. CZ-RAJ is situated in the Drahany Highland at a lower altitude and has a moderately warmer and drier climate than CZ-BK1. The reference evapotranspiration amount at each forest stand was computed from in situ data and was additionally used for quantifying the atmospheric evaporative demand from 2014 to 2016 (Table 1).

Eddy Covariance and Ancillary Measurements
At each station, the EC system consisted of an infrared gas analyser (LI-COR, Lincoln, NE, U.S.A.) and an ultrasonic anemometer (Gill, Hampshire, U.K.) measuring at a 20 Hz frequency (models are specified in Table 2). Each system was mounted on a tower several meters above the forest canopy. The EC measurements were complemented by an extensive set of sensors to collect the required auxiliary meteorological data. Measurements included the air temperature (Tair) and humidity profile with EMS33 temperature and humidity sensors (EMS, Brno, Czech Republic) and precipitation using Precipitation Gauge 386C (Met One Instruments, Grants Pass, OR, U.S.A.). Measurements for the incoming photosynthetic active radiation was made using LI-190R Quantum Sensor (LI-COR, Lincoln, NE, U.S.A.) at Bílý Kȓíz and EMS12 sensor (EMS, Brno, Czech Republic) at Rájec. Additionally, profiles of soil moisture were measured by using the CS616 (Campbell Scientific, Inc., Logan, UT, U.S.A.) sensors. An overall description of the instrumentation at the study sites is given in Table 2.
The VPD for each site was computed from Tair and relative humidity (RH) as follows [39]: where SVP (Pa) is the saturated vapour pressure given as follows: Tair 237.3+Tair (2)  At both study sites, the EC data from the main growing season (May-September) of 2014-2016 were used for the analysis. Gas analyser measurements of the density of the water vapour in the air and ultrasonic anemometer measurements of three-dimensional wind components and sonic temperature were made [35,36]. The eddy covariance processing was performed, using an open-source software, EddyPro 7.0.6 (Li-COR, Lincoln, NE, U.S.A.). The most recent methods for flux corrections, conversions, and a thorough quality control scheme as proposed by [40] were also applied. This process involves despiking and raw data statistical screening, basic quality checking (QC) of turbulent fluxes (flux stationarity and integral turbulence characteristics tests), coordinate rotation using the planar fit method [41], spectral correction [42][43][44], detecting and compensating for time lags of signals from the ultrasonic anemometer and the gas analyser, footprint estimations and calculating half-hourly fluxes. In EC data post-processing, determining periods with low turbulent mixing is a critical step [35]. Flux measurements over periods with insufficiently developed turbulence, i.e., low friction velocity (u * ) need to be detected and filtered out. This filtering procedure assured the exclusion of CO 2 measurements not representative of the biotic flux, i.e., net ecosystem exchange (NEE) [31,45,46]. The R software [47] package 'REddyProc' [48] was used to gap-fill the EC data at both sites, using marginal distribution sampling (MDS) [46]. This way, all missing values were replaced by the average value estimated under similar meteorological conditions within a specific time window. If no similar conditions were available within the starting time window of about 7 days, the window was increased. NEE was partitioned into GPP and ecosystem respiration (R eco ). The flux partitioning approach by [49] using daytime data was applied to estimate half-hourly GPP (µmol m −2 s −1 ) values. The half-hourly GPP values were aggregated to obtain daily and monthly GPP sums.

Soil Water Content Simulations
Since the soil volumetric water content (SVWC) was not measured throughout the entire period of our analyses (the measurements started at the beginning of 2016), we used a simulated daily SVWC instead. We applied the soil water balance model R-4ET, an R package for Empirical Estimate of Ecosystem EvapoTranspiration as used in [50]. The calibration of the soil water balance model was carried out using Bayesian statistics implemented in the R package BayesianTools [51] with a Differential-Evolution Markov chain Monte Carlo sampler. In Bayesian calibration, the model input parameters are iteratively updated to provide a probability distribution of the calibrated parameters that represents the uncertainty in the measured data and in the model structure. The simulations were repeated 4.8 × 10 6 times, where the first 1.8 × 10 6 runs were based on the prior distribution, while the remaining 3 × 10 6 were constrained by the posterior distribution resulting from the first set of runs. A number of iterations treated as burn-in were set to 0, and the thinning parameter determining the interval in which values are recorded was set to 1. The high number of iterations was necessary to attain a good input parameter convergence with a narrow distribution. For inspecting the convergence, we used Gelman diagnostics with a criterion for the potential scale reduction factor being less than 1.2 for all parameters [52,53]. The final selection of parameters from their probability distributions was conducted by using a maximum a posteriori probability estimate, i.e., the value representing the mode of the posterior distribution.
The main model input variables included meteorological data and the leaf area index (see [50], for more details). The soil was stratified into 12 layers, which increased in thickness with increasing depth down to 3 m. Note that the soil layers were selected in the way that the layers used for optimization matched the depths of all sensors with an extra ±2.5 cm, considering the volume measured by the sensors. At both sites, the available SVWC measurements were available throughout the main part of the root zone. At CZ-BK1, the sensors CS616 (Campbell Scientific, Inc., Logan, UT, U.S.A.) were placed at 0.05, 0.1, 0.22, 0.34, and 0.42 m soil depths and measured since 2016. At CZ-RAJ, the same type of sensors were also placed at 0.05, 0.1, 0.2, 0.5, and 0.8 m depths and measured since 2016. The soil parameters, including SVWC at saturation, field capacity, wilting point, and saturated hydraulic conductivity, were optimized at all 12 depths [50]. Additional single parameters relevant for the SVWC simulations that were optimized included the following: rooting depth with the Beta parameter describing the root profile shape [54], surface resistance and the degree of isohydricity [55], the water interception capacity of leaf and bark area, and curve number representing a runoff parameter [50]. The uniform distribution of priors was used, where the lower and upper limits were set to be within ±50% of the values based on field measurements of the wilting point, field capacity, and saturated water content and ±100% for the remaining parameters that were estimated from the literature or from previous anecdotal analyses. To provide the model with sufficient spin up time to stabilize and ensure reliable and robust parametrization, the simulation was initiated at the beginning of 2010 with initial conditions of SVWC set to the field capacity estimated from soil texture (note that 2010 was one of the wettest years at both sites, according to a computed standardized precipitation-evapotranspiration index over the region). The overall simulation was done for the period 2010-2019 where the observed data for Bayesian calibration spanned the period 2016-2019.
The simulated SVWC averaged over all depths yielded a root mean square error of 0.037 m 3 m −3 and 0.017 m 3 m −3 for CZ-BK1 and CZ-RAJ, respectively, suggesting realistic SVWC estimates.

Drought Stress Determination
According to [56], drought refers to conditions characterized by low available soil moisture and high atmospheric VPD. As such, in quantifying the dry periods in our study, we used the Standardised Precipitation-Evapotranspiration Index (SPEI), which takes into account both precipitation and potential evapotranspiration for the determination of drought over the main growing season. The SPEI was calculated for various lags (1, 3, 6, 12, and 24 months) from monthly records. We used the period 1981-2010 as the baseline period for the SPEI computation. The computation of SPEI was performed according to [57,58], using the R package SPEI [59]. The SPEI represents an anomaly in the climatological water balance, which is given by the difference between precipitation and reference evapotranspiration. The input data for the SPEI derivation were obtained from climate reanalyses ERA5 monthly averages available at 0.25 • spatial resolution, where the reference evapotranspiration was computed from air temperature at 2 m, air dew point temperature in 2 m, wind speed in 2 m converted from the original 10 m height using a logarithmic profile law, and shortwave incoming radiation, following the methodology by [60]. In determining the occurrence of dry and wet events in our study, SPEI classification based on [58] was used (Tables 3 and 4). Light response curves (LRC) of daytime GPP were fitted at half-hourly time resolution using the logistic sigmoid approach by [61]: where PAR is the photosynthetically active radiation, α is the apparent quantum yield in mol (CO 2 ) mol −1 (phot.) that describes the effectivity of photosynthesis at low light conditions, and the GPP max is the asymptotic maximum assimilation rate at the light saturation point in µmol m 2 s −1 . The fitting was done separately for half-hourly values from years with normal conditions (2014 and 2016) and the year with drought stress conditions (2015). To eliminate night-time measurements, we used a PAR threshold 10 µ, i.e., light compensation irradiance of 10 µmol m 2 s −1 , representing light compensation point of NEE. Since VPD and soil moisture affect GPP [49,62], the effects of VPD and SVWC on the LRC residuals (GPP values after the removal of the strongest PAR dependence) were analysed for both spruce forest sites within the period under study. A piecewise regression was performed on the LRC residuals versus VPD and SVWC to determine the site-specific environmental stress thresholds at which GPP declined during the growing season periods of normal and drought affected years.

Piecewise Regression Analyses for the Assessment of Drought Effect
Since the LRC model does not account for changes in both VPD and SVWC, piecewise regression of residuals against these environmental variables were analysed. This allowed to determine the response of the studied spruce forest ecosystems with contrasting climates to severe drought conditions over the growing season for normal and drought-affected years. A piecewise regression analysis as part of the R package in 'segmented' [63] and the Davies test [64] were used to detect the breakpoints in the regression and also to test for the significant differences in slope parameters of a plot of residuals from the LRC model to drought conditions (high VPD and low SVWC).

Differences in Meteorological Conditions at the Experimental Sites
During May-September from 2014 to 2016, the year 2015 was exceptionally dry compared to 2014 and 2016 (Table 4). Additionally, 2015 was characterized by significantly higher mean VPD values (p < 0.001; Table 5) during the main growing season as compared to the same period in 2014 and 2016 across both sites. Hence, the mean VPD and SVWC values were the best meteorological indicators of drought, while Tair was comparable for years 2015 and 2016. Consequently, the high atmospheric evaporative demand (high VPD values) and low soil moisture (low SVWC values) recorded in 2015 across both sites during the study period further indicated the dryness stress experienced in the main growing season of 2015. There were observed statistically significant differences (p < 0.05; Mann-Whitney-Wilcoxon ranksum test) in the VPD and SVWC values between the normal and drought regimes across both sites (Figure 1). In addition, VPD was higher and SVWC lower at CZ-RAJ than at CZ-BK1, especially during the drought-affected year (DY).

Light Response Curves of GPP at Different Climates
The parameters of LRC for normal years (NY) 2014 and 2016 were compared with those for DY for both experimental forest sites (Figure 2; Table 6). Generally, the apparent quantum yield (α) at CZ-BK1 was found to be higher than that in CZ-RAJ. During the years with normal conditions, α in CZ-BK1 was observed to be 12% higher than in CZ-RAJ when PAR was less than 500 µmol m −2 s −1 . Moreover, during the DY of 2015, α at CZ-RAJ further declined by 25% as compared to that in CZ-BK1.
The maximum gross primary production at light saturation (GPP max ) for CZ-BK1 was found to be higher than that for CZ-RAJ during the entire study period. During the NY period, GPP max in CZ-BK1 was 34% higher than in CZ-RAJ. However, during the dry year period, there were significant reductions in the GPP max values recorded at both spruce forest stations with (18% decline in CZ-BK1 and 17% decline in CZ-RAJ) as compared to those recorded for the years with normal conditions. Additionally, it was further observed that the GPP max value recorded at CZ-RAJ during the DY period was 33% lower than that in CZ-BK1. Interestingly, the GPP max value recorded in CZ-BK1 during the DY period was still higher than the GPP max value recorded at CZ-RAJ during the years with normal conditions.  Table 6. Light response curve parameters for the normal (2014 and 2016) and dry (2015) years for the wet and dry climates (CZ-BK1 and CZ-RAJ respectively) within the main growing season period of May-September. The apparent quantum yield (α), the maximum gross primary production at light saturation (GPP max ) and the coefficient of determination (R 2 ) are also shown.

Response of Light Response Curve Residuals to VPD and SVWC at Different Climates
The applied LRC model reflects only the light response of GPP and thus, its residuals can be used to assess the influence of additional meteorological parameters (i.e., VPD and SVWC) that are known to affect GPP (Figures 3 and 4). The residual analysis confirmed that the residuals consistently tended to be more negative with increasing VPD and decreasing SVWC (dry conditions). VPD had a significant and stronger effect on GPP than SVWC across both sites, as seen from the piecewise regression analysis, using the residuals from the LRC for the years under study. However, during the DY period, both VPD and SVWC had significant effects on GPP at CZ-RAJ. For CZ-BK1, there was a minimal effect of SVWC on GPP during the years with normal conditions. Generally, the relationship between the residuals and changes in VPD and SVWC revealed a biphasic response to drought, except for the DY period in CZ-RAJ in Figure 3. All the breakpoints from the piecewise regression were found to be statistically significant (Table 7). However, steeper slopes after the VPD breakpoint values from the initial slope were mostly observed for both DY and NY periods in CZ-RAJ and only for DY period in CZ-BK1. This shows that during all the years under study, GPP decreased at a faster rate after the breakpoint in CZ-RAJ than in CZ-BK1, except for the DY period when there was a significant decline in GPP at a faster rate after the breakpoint across both sites.

Impact of Drought on Ecosystem Carbon Fluxes at Monthly Timescale
The monthly averages of daily sums of GPP and R eco were analysed separately to assess whether the observed net ecosystem productivity (NEP) reduction during drought was mainly a consequence of increased R eco , reduced GPP, or contribution of both across the forest ecosystems (Figures 5-7).      It was observed that there was a significant decline in the total GPP by 14% and 6% during the main growing season period of the dry year for CZ-RAJ and CZ-BK1, respectively. There were also observed significant statistical differences (p < 0.001) in the mean monthly GPP values recorded for July and August during the dry year as compared to the years characterized by normal conditions at both forest sites, especially for CZ-RAJ ( Figure 5). Moreover, this significant decline in the mean monthly GPP values at both spruce forest stands during the months of July and August coincided with high daily mean VPD values (>12 hPa) and low SVWC values (<0.16 m 3 m −3 ) as compared to the same period for the adjacent years with normal conditions. However, an increase in the monthly mean R eco values was observed for July to September of the dry year as compared to the adjacent normal years within CZ-BK1 but with no significant statistical differences, except for the month of August ( Figure 6). These three months within 2015 were also characterized by high mean Tair (>20 • C) and low mean SVWC (<0.19 m 3 m −3 ) as compared to the two other adjacent years. In addition, during the DY period within CZ-RAJ, the monthly mean R eco values were found to have declined significantly, as compared to the adjacent NY periods. We also observed a 38% significant decline in the mean monthly NEP values during the dry year period for the dry spruce forest in Rájec, compared to a 12% decrease in the mean monthly NEP values for the humid spruce forest in Bílý Kříž (Figure 7). There were also statistically significant differences (p < 0.001) in the mean monthly NEP for the months of July and August across both spruce forest stations. During these months, the observed mean monthly NEP values largely declined during the DY period, as compared to the same period during the NY, especially in CZ-RAJ.
To summarize, the large decline in GPP, R eco and NEP during the dry year period (especially from July to August) for CZ-RAJ showed that the impact of the drought was more severe in CZ-RAJ than in CZ-BK1. However, during July-August of the DY period, R eco at CZ-BK1 significantly increased as compared to that in NY.

Discussion
The study sought to assess the different effects of the 2015 summer drought on wet and dry spruce forest ecosystems. Our findings corroborate the earlier published results by [27][28][29] that the year 2015 was characterized by severe drought conditions across Europe ( Figure 1 and Table 4). Furthermore, it was found that the months of July-August of 2015 were the most affected by high mean Tair (>20 • C), high mean VPD (>10 hPa) and low mean SVWC values (<0.19 m 3 m −3 ) across both spruce forest stands. This shows that the drought effect in 2015 was severe in the months of July-August. Moreover, the results of this study also show that the rate of the forest ecosystem photosynthesis was significantly reduced (likely through the immediate closure of the stomata to protect the tree from desiccation) by the high VPD, Tair and soil water deficit during these two months in 2015. As such, there was an observed strong decline in both forest GPP and NEP across both spruce forest stands during the DY period of 2015, especially in CZ-RAJ compared to CZ-BK1 (M). This strong decline in forest GPP and NEP at CZ-RAJ was mainly due to the high atmospheric evaporative demand coupled with the low SVWC conditions experienced at this forest ecosystem, as compared to that in CZ-BK1 [30]. However, due to the humid climatic conditions at CZ-BK1, an increase in Tair during the months of July-August in 2015 only aided the increase in the kinetics of the enzymes participating in microbial decomposition and root respiration under warm conditions, thereby increasing the overall forest ecosystem respiration [65][66][67][68][69][70][71].
Additionally, results from the piecewise regression analysis using the residuals from the LRC model highlight the effects of both VPD and SVWC on forest ecosystem GPP. This is because both photosynthesis and transpiration are mediated by stomatal conductance, which are affected by these environmental variables. However, a steeper decline in forest ecosystem GPP (Figure 3 and Table 7) was observed with high VPD values across both forest stands, even under non-drought years (when SVWC was non-limiting). This is consistent with recent studies, which show that high VPD values aggravate drought effects in forests, due to the abrupt changes over very short timescales within a day, even without dry soil conditions [7,72,73]. This further explains the strong suppression of GPP by high VPD on even non-drought years (2014 and 2016), as the SPEI (for determining dryness) only captures changes in the contributing environmental drivers (temperature and soil moisture) for longer time scales of weeks or months [74,75]. Thus, we would recommend to include the influence of these abrupt changes in VPD on the rate of photosynthesis with SVWC limitations in future LRC models when analysing the impact of drought on GPP, especially for different forest ecosystems that are exposed to regular strong edaphic droughts [49,62,76,77].

Conclusions
The study shows that atmospheric constraints increase the vulnerability of the Norway spruce forest to the severity of drought, especially at sites with a moderately dry climate that are characterized by precipitation that is typically equal or smaller than the atmospheric evaporative demand. It also further highlighted the strong influence of VPD on carbon uptake, which was further worsened by the decline in soil moisture. The effect of SVWC on GPP was especially noticeable during severe drought conditions within the DY period. Consequently, with regards to climate change, our results suggest that elevated temperatures will further exacerbate the drought impacts on forest (Norway spruce) ecosystems at sites with precipitation levels equal or smaller than the atmospheric evaporative demand, such as CZ-RAJ. The decline in GPP and NEP in 2015 found in our study thus questions not only the sustainable productivity, but also the existence of Norway spruce per se in such areas, considering the prolonged period of drought in future climatic conditions. The results of this study may help decision makers to quantitatively assess the performance of Norway spruce in future climatic conditions.

Data Availability Statement:
The data presented in this study are available on FLUXNET and could also be made available on request from the corresponding author.

Acknowledgments:
The authors wish to thank Radek Czerný who assisted in the eddy covariance data processing.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.