Changing climate sensitivity of secondary growth following extreme drought events in forest ecosystems: a global analysis

Understanding tree-response to extreme drought events is imperative for maintaining forest ecosystem services under climate change. While tree-ring derived secondary growth measurements are often used to estimate direct and lagging drought impacts, so-called drought legacies, underlying physiological responses remain difficult to constrain across species and site conditions. As extreme droughts may alter the functioning of plants in terms of resource allocation being shifted towards repair and physiological adjustments, climate control on growth may consequently be altered until physiological recovery is completed. In this context, we here advance the concept of drought legacy effects by quantifying ‘functional legacies’ as climate sensitivity deviations (CSD) of secondary growth after droughts, i.e. temporary alterations of climate-growth relations. We quantified climate sensitivity deviations after extreme drought events by applying linear mixed-effects models to a global-scale, multi-species tree-ring dataset and differentiated responses by clades, site aridity and hydraulic safety margins (HSMs). We found that while direct secondary growth legacies were common across these groups, responses in post-drought climate sensitivity deviations were nuanced. Gymnosperms showed weaker coupling between secondary growth and the dominant climatic driver after droughts, a response that was narrowed down to gymnosperms with a small HSM, i.e. risky hydraulic strategy. In comparison, angiosperms instead showed stronger coupling between secondary growth and the dominant climatic driver following droughts, which was narrowed down to the angiosperms growing in arid sites. These results are consistent with current understanding of physiological impairment and carbon reallocation mechanisms, and the distinct functional responses suggest functional legacies quantified by climate sensitivity deviations is a promising avenue for detecting and thus studying physiological mechanisms underlying drought-responses in tree growth on large scales.

Trees exposed to droughts need time to recover their physiological functionality [17,18], resulting in legacy periods that can last from a few days [19] to multiple years [20]. Impaired physiological function under severe drought results in reduced gross primary productivity as pre-drought carbon assimilation rates cannot be maintained [18,20,21]. Since respiratory costs rise strongly with increasing temperature, net primary productivity (NPP) is even more reduced during droughts [18], especially as secondary growth [16]. In this context, tree-ring widths serve as a valuable proxy for tree growth performance in relation to droughts since NPP is linked to secondary growth [22]. Indeed, legacy periods of reduced growth identified from tree-ring measurements have been found to be widespread [13,20,[23][24][25] and lasting 1-4 years after droughts [6,13,26,27].
Two of the most important and frequently discussed physiological mechanisms involved in trees' responses to drought are hydraulic failure and carbon starvation [30,33,34,[53][54][55]. Hydraulic failure is the partial or complete loss of water conductivity due to conduit embolism and is more likely in anisohydric species which maintain high stomatal conductance despite increasing water deficit [56]. Carbon starvation is caused by depletion of nonstructural carbohydrates (NSCs) and is more likely in isohydric species which in contrast to anisohydric species downregulate stomatal conductance to conserve water with the consequence of reduced carbon assimilation [7]. Because of the need to balance risks of hydraulic failure and carbon starvation, hydraulic strategies have evolved with a strong climatic adaptation [30,32,33]. The hydraulic safety margin (HSM) is an informative metric indicating how risky species' stomatal conductance strategies are, and is calculated as the difference between the daily minimum water potential under non-extreme conditions and the water potential deficit required for a 50% loss of conductivity (ψ 50 ) [33,57,58]. Large HSM-values are related to isohydric strategies and imply a low risk of experiencing cavitation-inducing conduit tensions, while small HSM-values are related to anisohydric strategies and imply that cavitationinducing conduit tensions in the xylem are frequent [30,33,55,59].
Extreme droughts necessitate recovery of damaged hydraulic functionality or altered carbon allocation as (re)growth of foliage, roots, or both are prioritized over secondary growth [60]. As a consequence, the relationship between secondary growth and climate can be temporarily altered [61]. This temporal change in climate sensitivity, climate sensitivity deviation (CSD), is driven by physiological mechanisms induced by drought rather than by nonstationary climatic conditions, which occur naturally over longer temporal scales [62,63]. Compared to the primary variable of interest (tree-ring widths as proxy for NPP), CSDs are identified at the level of functional relations of radial growth with climate variables, with strengthened or weakened relationships implying more or less climate control on tree-ring widths [61]. CSDs are therefore distinctly different from the traditional legacy effect based on direct secondary growth which has been extensively reported in recent years [13,20,[23][24][25].
A weakening in post-drought climate-growth relations (decoupling between climatic variables and radial growth) would be expected based on the observation of preferential carbon allocation to other metabolic costs or non-secondary growth as well as due to embolism [60,61], with greater or lesser intensity of response depending on abiotic and genotypic differences. Indeed, such changes in climate-growth relationships after droughts have been reported on local scales [44,51]. However, a global perspective across regions and taxa is currently lacking. Complementing the recent surge in scientific interest on legacy effects in direct secondary growth, additional analyses for changes in climate sensitivity (CSD) can identify more nuanced response patterns. As CSDs hint at underlying mechanisms not detected in direct secondary growth responses, they may result in a more comprehensive understanding of trees' response to droughts.
To this end we here selected sites from a global tree-ring archive, the International Tree-Ring Databank (ITRDB), which are sensitive to droughtindices and quantify 'functional legacies' as lagged CSDs after extreme droughts. In conjunction with gridded climate products, the ITRDB allows for comparison between clades and levels of site aridity, while available data of HSMs for many of the analyzed species also allows comparison between hydraulic strategies. Thereby, factors that convey disadvantages in the face of extreme droughts can be identified. We thus address the following hypotheses: (a) extreme drought events that elicit ecological responses are succeeded by weakened (decoupled) drought-index sensitivity of secondary growth (CSDs), and (b) functional legacies detected by post-drought CSDs vary depending on the factors influencing drought response: clade, site aridity and HSM.

Tree-ring data
We used the harmonized version of the ITRDB presented in Zhao et al [64] containing 4207 sites, available at the National Climatic Data Centre (www.ncdc.noaa.gov/paleo-search/study/25570). Although Europe and North America are overrepresented, the dataset contains tree-ring data from all continents but Antarctica. The raw ring-width data was post-processed and quality-checked following standard protocols (see SI Methods for details). Furthermore, we limited chronologies to those: (a) with a length of at least 25 years after 1950, to match the temporal range of climate data; (b) whose meta-data specified species or genus; and (c) with plausible site coordinates which reduced the 4207 sites to 2805 (SI figure 1).

Climate data
Monthly gridded mean, minimum and maximum temperature data since 1950 at a spatial resolution of 0.5 • were obtained from the Climate Research Unit (CRU) TS 4.04 dataset [65] and equally resolved precipitation data were obtained from the Global Precipitation Climatology Centre (GPCC), version 2018 [66]. Precipitation data from GPCC was used over CRU due to better capacity to predict growth of the ITRDB sites (SI figure 2), likely because of the higher station density [66]. We only used climate data after 1950 because of low station density prior to that year. To conform to local topographical effects of the tree-ring sites, we statistically downscaled these gridded products from 0.5 • to 30 arc seconds resolution [67] using CHELSA version 1.2 [68]. We calculated potential evapotranspiration (PET) [69], from which we derived climatic water balance (CWB) as CWB = P−PET [70] as well as integrations of the standardized precipitation-evapotranspiration index (SPEI) over 4-12 months (SPEI4-SPEI12, from here on collectively called SPEIx, SI figure 3) [71]. We avoided shorter SPEI time-scales as tree-rings typically feature longer response times to drought indices [72,73] and lower time-scales have a higher probability of reaching extremes due to frequent variation [74].

Data processing
A significant positive correlation between ring-width indices (RWIs) and drought index was required for sites to be included to avoid spurious results [13]. Site-dependent determination of SPEI time-scale was preferred over a generalized single combination to account for site-dependent growing seasons and thus SPEI time-scale responses [45,72,73,75]. This method of selecting climate dependencies for secondary growth led to considerably more sites and drought events being retained compared to if a generalized combination across all sites was used (SI figure 5). Hence, the month and SPEI time-scale with the highest significant correlation were selected for each site (figure 1, see SI Methods for details, SI figure 3), resulting in 1236 sites being excluded and 1569 sites being retained for further analyses.
We applied a novel approach for defining severe droughts by necessitating that climatological droughts have corresponding negative growth anomalies, i.e. extraordinarily low RWI, to identify droughts that led to growth decline. Not every drought will lead to growth decline (SI figure 4 [76]) due to differences in drought timing and intensity [73,77], nor do all growth declines correspond to droughts. This is particularly true when using a probability-based drought index such as SPEI which, Table 1. Number of unique sites and species per top-and bottom-level group. Total number of unique sites and species are 793 and 108, respectively. Groups beginning with asterisk ( * ) are hydraulic safety margin (HSM) related angiosperm groups which were not used in the final analyses due to being insufficiently represented. when negative, implies lower than average water availability, not necessarily water deficit [78]. Therefore, we defined climatological droughts as the cooccurrence of a −2 standard deviation in SPEI [79] and negative CWB for the month and SPEI time-scale selected for each site, while for the corresponding ecological responses, we scanned all chronologies for negative anomalies using a recently published pointer year detection method, the bias-adjusted standardized growth change method (BSGC) [80]. The BSGC method relies on probability density functions of annual growth changes and defines abnormal growth changes outside the 95% confidence interval. This method removes the arbitrary thresholds used by traditional methods to identify pointer years, while additionally resulting in a more precise pointer-year detection rate for both existing and artificially generated tree-ring data [80,81]. The combination of BSGC-determined negative pointer years and climatological droughts constrained our analysis to all droughts that featured extraordinary growth changes, ensuring significantly lower RWI-values compared to only using the climatological thresholds (Wilcoxon rank-sum test, p < 0.001, SI figure 4). Despite this comparatively strict threshold, at least one drought event was found in 793 sites out of the remaining 1569, while none were detected in the other 776 which were therefore not considered further. In summary, 793 sites covering 108 unique species (see SI table 1 for a detailed list): (a) passed statistical quality and representativeness; (b) had at least one significant positive SPEIx-month correlation; and (c) possessed at least one extreme drought event and post-drought year after 1950, and were thus retained for further analyses (figure 1). We used the clade, site aridity and HSM to divide the data into groups (table 1). Clades were assigned based on genera, and sites were categorized as arid or wet based on whether the mean annual CWB (1950-2015) was negative or positive, respectively. The separation between arid and wet sites was validated through a sensitivity analysis where all sites with mean annual CWB of ±100 mm and ±250 mm were removed (see SI Methods, SI figure 10). ψ 50 HSM data were taken from Choat et al [33] and the TRY plant trait database [82], and covered 30 of our species (28%), corresponding to 514 of the utilized sites (65%). We used the static HSM-value of 1 MPa as the difference between a small (<1) and large (>1) HSMs as suggested in previous studies [33,83]. However, the cover of HSM-values for angiosperms were insufficient to compare small and large HSMs (four vs one species with large vs small, respectively, table 1). Therefore, only gymnosperms were compared in regards to HSM, for which the data covered 457 sites (71%) and 25 species (30%). Because of the incomplete coverage of HSM-values and uneven distribution of sites for the species in the ITRDB, a sensitivity analysis for the 1 MPa border similar to that of mean annual CWB of site aridity, was not possible. The categorization resulted in four 'top-level' groups, based on clade and site aridity (SI figure 6), and six 'bottom-level' groups, where clade, site aridity and, for gymnosperms, HSM were crossed (e.g. gymnosperms in arid sites) (figure 2, table 1).

Statistical analyses
Based on previous studies, we set the analyzed postdrought (and thus maximum potential legacy) period to four years [13,26]. If subsequent droughts were detected within a post-drought period, the year of the new drought was considered as another drought year and the count of the post-drought period restarted (thus cutting the previous post-drought period short). To determine post-drought CSDs we used linear mixed-effect models (LMM ) with RWI as a function of SPEI indices and year-type. This model structure allows for determination of changes in direct RWI and sensitivity of RWI to SPEI in each postdrought year compared to all non-drought or nonpost-drought years, the last of which we from here on call 'other years' . Sites were assigned random intercepts to account for site-specific differences [84]. The LMM was applied on all data groups individually (SI figures 7 and 8), allowing comparison of postdrought effects between groups, thus corresponding to our second hypothesis. Further details on the LMM structure are presented in the SI methods. The analyses were done in R v4.2 [85] and RStudio v1.4.1743 [86] extended with the packages tidyverse [87], SPEI [88], raster [89], dplR [90], broom [91], nlme [92], patchwork [93], dendRolAB [94], sf [95], rnaturalearth and rnaturalearthdata [96].

Results
Linear mixed models (LMMs) applied on sites grouped by clade revealed significant and cladespecific changes in drought-index sensitivity during years following drought events. While some gymnosperms showed weakened drought-index sensitivity (negative CSD) after drought (figure 3(a) rows 5-9), some angiosperms showed strengthened drought-index sensitivity (positive CSD) (figure 3(a) rows 2-4). Consequently, across 839 drought events in 793 sites, divergent changes in drought-index sensitivity were found. Although post-drought RWI deviations also displayed diverging patterns, the majority of deviations for both clades were expressed as reductions ( figure 3(b)), which altogether resulted in a consistent RWI legacy lasting for three years (figure 3(b) row 1, SI table 2).

Gymnosperms
The largest group, gymnosperms, showed significantly negative CSD at two years after drought (figure 3(a) row 5, SI table 2), indicating weakened drought-index sensitivity of RWI to SPEI. In comparison, persistent RWI reductions lasted for three years after drought for the same group (figure 3(b) row 5, SI table 2). Differentiating gymnosperm sites by aridity did not reveal any significant CSDs, while the three-year RWI reduction remained for arid sites (figure 3(b) row 6, SI table 2) whereas it lasted for two years for wet sites (figure 3(b) row 7, SI table 2). The subset of gymnosperms which featured HSM values showed that species with small HSM displayed significant negative CSDs that persisted for the first two years after a drought (figure 3(a) row 8, SI table 2). The corresponding RWI legacy lasted for two years also for this group (figure 3(b) row 8, SI table 2). In comparison, no CSDs nor RWI legacy were observed for gymnosperms with large HSM. However, RWI deviated negatively at three years after droughts (figure 3(b) row 9, SI table 2).

Angiosperms
As a whole group, angiosperms did not display any significant CSDs after drought (figure 3(a) row 2, SI table 2), while the corresponding deviations in RWI showed a shorter legacy compared to gymnosperms at only one year after drought (figure 3(b) row 2, SI table 2). Differentiating angiosperms by site aridity revealed that those in arid sites showed positive CSDs for two years after drought (figure 3(a) row 3, SI table 2), indicating strengthened drought-index sensitivity of RWI, rather than weakened as detected for gymnosperms. The corresponding RWI legacy of angiosperms in arid sites showed an initial one-year reduction in RWI followed by a persistent overshoot for the subsequent three post-drought years (figure 3(b) row 3, SI table 2). In comparison, no significant changes in drought-index sensitivity were seen for angiosperms in wet sites, although RWI was reduced for one year in this group as well ( figure 3(b)

Discussion
Our analyses reveal that altered climate-growth relationships-climate sensitivity deviations (CSDs)-after extreme drought events prevail in angiosperms in arid sites and gymnosperms with small HSM (figure 3(a) rows 2 and 5), which together make up 37% of all sites (294 out of 793) across the spatial extent of the utilized ITRDB subsample. The weakened climate sensitivity of gymnosperms with small HSM supports our initial hypothesis of a functional post-drought legacy as it vanishes, i.e. returns to normal, two years later [61]. Although our results cannot conclusively prove what mechanisms are behind the CSDs, the observed positive CSDs of angiosperms in arid sites is consistent with NSCmobilization towards xylogenesis (see section 4.2.1), and therefore also fit within the homeostatic sensitivity framework [61].

Small HSM
Among gymnosperms, a significant negative CSD legacy was found for species that exhibited small HSMs, with a simultaneous two-year RWI reduction ( figure 3, row 8). These results are consistent with physiological impairment and consequently less efficient exploitation of precipitation after drought [18,97]. The overall vulnerability of species with small HSM to droughts is also consistent with other studies [13,32,53,[98][99][100], especially for gymnosperms which compared to angiosperms have less NSCs available in the xylem which are necessary for removing nearby embolism [101]. It has been suggested that trees with small HSM rely on other architectural and dynamic traits, such as rooting depth, sapwood capacitance and cell wall resistance to counteract hydraulic dysfunction in course of droughts [30,[102][103][104]. However, our results suggest these mechanisms to be insufficient for mitigating ongoing impacts of extreme drought in comparison to species with large HSM.

Large HSM
Compared to gymnosperms with small HSM, those with large HSM did not display any significant CSDs nor immediate RWI legacy (figure 3, row 9), suggesting that regulation of stomatal conductance to maintain a HSM of >1 MPa may generally be sufficient to avoid hydraulic failure after extreme drought events [53,102]. However, this cutoff should be taken with caution as there are problems with HSM values originating from measurements of daily minimum water potentials being higher (less negative) than in reality [105]. The consequence is that species may be considered to have smaller HSM than they have in reality, which then implies the 1 MPa cutoff used in our analyses may in fact be higher [105]. However, irrespective of this issue, the diverging responses seen for gymnosperms with small and large HSM remains. The delayed RWI reduction at three years after drought may reflect acclimation via physiological downregulation [60,106,107], as suggested by the concurrent maintained climate-growth relations. Carbon allocation towards rebuilding NSCstores works on a shorter time-scale and is unlikely to be involved [97,108].

Site aridity
While gymnosperms with small HSM stand out among gymnosperms with a consistent CSD legacy, gymnosperms in both arid and wet sites also show multi-year RWI reductions ( figure 3(b), rows 5-8), thus agreeing with many other studies that multiyear RWI legacies are widespread [13,20,[23][24][25][26]109]. Furthermore, the gymnosperms in wet sites showed larger reductions in RWI compared to those in arid sites, which is in line with other studies [8,36,38,110] and can be explained by resource economics. Plants adapt growth rates depending on average growing conditions [111], resulting in strategies with faster growth on wet sites [33,112,113] and associated higher drought susceptibility [17,20,29,40,100,102] via larger crowns and lower cavitation resistance due to wider xylem conduits and lower wood density [36,114]. Consequently, reductions in radial growth may be larger compared to sites with lower water availability, while the inherent higher average water availability also appears to convey a faster growth rate recovery. The lack of CSD for gymnosperms in wet sites is surprising as the above reasons for larger impacts on growth also suggests a more vulnerable physiology. A possible explanation for this may be that site wetness convey a lesser influence on the physiological vulnerability of the trees compared to their stomatal regulation, wherefore CSDs are detected when HSMs are compared. Conversely, and following the same principles of resource economics, gymnosperms in arid sites avoid physiological impairment by an evolutionarily adapted physiology, such as more conservative growth rates [38,111,115], bimodal xylogenesis [116], better xylem structural features such as higher cell wall resistance and wood density [29][30][31], and a more developed root system [48,101,117,118]. While RWI legacies still occur, they appear less severe compared to those in wet sites.

Potential conflation with masting investment
The CSD in the second year after drought for gymnosperms as a whole group is not considered consistent with a CSD legacy because of the discontinued response. However, the emergence at two years after drought may potentially indicate associated masting investments. If cued by warm conditions [119] and with negative impacts on secondary growth due to internal resource competition [120], masting investment can theoretically also lead to single-year CSDs. However, species differ in reproduction strategies, cues, growth sensitivity and maturation times, and consequently are assumed to be averaged out within the groups in this study. Additionally, we want to stress that sites with gymnosperms with small HSM constitute more than one third of the drought events recorded on sites with gymnosperms in this study, wherefore part of the observed CSD in the second year after drought might originate from those. Although climatic cues and tree-ring-sensitivity of masting are unclear, a masting-related CSD may be seen for Ponderosa pine (Pinus ponderosa) which has a twoyear cone maturation time [121], and appears as a non-significant negative CSD in conjunction with a negative RWI deviation (SI figure 9).

Arid sites
The positive CSDs for angiosperms were contrary to our expectations as we hypothesized weakened sensitivity of RWI to climate following drought. Differentiation by site aridity showed that the positive CSD legacy was related to arid sites, which together with the successive RWI overshoot ( figure 3 row 3) indicate a physiological mechanism of enhanced carbon allocation towards xylogenesis after droughts. While removal of emboli from conduits through NSC mobilization is one mechanism for regaining hydraulic integrity [101,[122][123][124], especially in angiosperms [125], rebuilding xylem may be a more common and effective mechanism [7,17,18,97]. By xylogenesis, RWI would be assumed to follow the observed pattern of an initial impairments followed by overcompensation. The first post-drought year with low RWI is additionally strongly influenced by extraordinarily negative SPEI values (SI figure 7(c)), equaling a more limiting climate which by Liebig's law of the minimum strengthens the coupling between climate and RWI [126]. The subsequent RWI overcompensation could also be enhanced by decreasing competition due to drought-induced mortality of neighboring trees [52]. Such details of the sampling sites are unknown, and if mortality patterns play a role, they would contribute to a survival bias and consequently inflate the positive drought response. While droughts are generally connected to negative impacts on growth, overshoot growth after droughts has also been reported [13,21,27,127,128].

Wet sites
Angiosperms in wet sites, contrary to those in arid, do not show RWI overshoot after drought, but a slightly stronger RWI reduction in the first year. This would be expected given that trees in wet sites preferentially try to rebuild canopies [38,98] or expand root systems [129], either because of a lacking drought adaptation or because of a requirement to maintain current structure.

ITRDB
While we find indications of deviating climate sensitivities later than two years after droughts, the significant changes are limited to the first two ( figure 3(a)), indicating the physiology of most trees recovers within two years. Although this study takes a broad perspective on physiological effects that are highly dependent on micro-climates, the lack of longer lasting CSDs and persistent physiological impairment may be counterintuitive in the context of our current understanding of drought vulnerability of forest ecosystems. In this context, a possible complication is the well-known inherent survival bias of the ITRDB [64], especially for arid and warm climates [130], which is a natural consequence of sampling protocols mostly aimed at climate reconstruction [131] resulting in a sampling bias towards the most vital and dominant trees. While our results on one hand show CSDs are present even in potential survivors, this bias may also affect the observed duration of legacy effects. The increasing momentum of ecological studies involving tree-ring data [132] will hopefully lead to datasets of comparable size to that of the ITRDB in the near future which will counter biases like this.

Drought indices
Interpretation of the results are dependent on the choice of drought indices. Choosing site-dependent combinations of SPEI time-scales and months resulted in the climatic dependencies that best explain the respective sites' secondary-growth. The considerably lower number of sites and species showing significant positive correlations with fixed time-scales and seasons for SPEI (SI figure 5) support our preference for a flexible approach [72]. It is also reasonable to expect that droughts detected in the most relevant climate variable for inter-annual growth variability also leads to detecting the most relevant and impactful droughts. However, the consequence is a heterogeneous mixture of drought seasonalities (month) and antecedent climates (SPEI time-scale) (SI figure 3), both of which are vital if examining for related vulnerabilities and impacts [11,36,127], but which is not the aim of this study. Further, it is likely that using generalized SPEIx-month combinations leads to a selection bias which predominantly filters for sites from within the same bioclimatic regimes, potentially also explaining the difference in results in comparison to using flexible SPEI.

Conclusion
We have shown that functional legacy effects in the form of climate sensitivity deviations (CSDs) can serve as a valuable window to discover and differentiate underlying mechanisms of how trees' functionality responds to severe drought events. CSDs reveal that trees growing under different abiotic conditions and with different genotypic features can differ in drought response characteristics, irrespective of direct impacts on secondary growth. CSDs therefore further nuance impacts of droughts on trees compared to the traditionally solely reported change in secondary growth. Although CSDs alone can merely identify functional drought legacies, combining our results with current ecological understanding can expose relevant mechanisms. More specifically, the results of this study hint at diverging responses related to resource economics where trees that grow under drought-prone conditions are less affected, while trees growing under less drought-prone conditions are more affected and respond in line with maintaining a fast growth strategy. A strong difference is seen between clades, with angiosperms in drought-prone conditions indicating an adapted drought-response which could be due to increased focus on xylogenesis which translates into positive CSDs, while gymnosperms with small HSM may be vulnerable to drought events as the negative CSD suggests impaired physiological integrity. Negative growth legacies also differ between clades with it lasting 1 year for angiosperms compared to 2-3 years for gymnosperms.
Still, more detailed species-specific studies using additional data streams, especially from experimental ecophysiology, are vital to conclusively link CSDs to physiological mechanisms and find intraspecific differences. Likewise, confounding factors such as how masting events impact CSDs should be studied to avoid conflating signals, also in regards to changing sensitivities mediated by climate change (e.g. [133]). The different pathways between drought impact and response implies unequal drought effects, which in turn require different frameworks for evaluating fitness under climate change. In conclusion, the presented statistical quantification of CSDs across a wide range of ecological settings indicates new research avenues in the context of climate change impacts on tree performance.
All data that support the findings of this study are included within the article (and any supplementary files).
paper. C Z, A B and A R contributed input to study design and paper writing.

Conflict of interest
Authors declare no conflict of interest.