Adaptations of white spruce to climate: strong intraspecific differences in cold hardiness linked to survival

Abstract Understanding local adaptation of tree populations to climate allows the development of assisted migration guidelines as a tool for forest managers to address climate change. Here, we study the relationship among climate, a wide range of physiological traits, and field performance of selected white spruce provenances originating from throughout the species range. Tree height, survival, cold hardiness, hydraulic, and wood anatomical traits were measured in a 32‐year‐old common garden trial, located in the center of the species range. Provenance performance included all combinations of high versus low survival and growth, with the most prevalent population differentiation for adaptive traits observed in cold hardiness. Cold hardiness showed a strong association with survival and was associated with cold winter temperatures at the site of seed origin. Tree height was mostly explained by the length of the growing season at the origin of the seed source. Although population differentiation was generally weak in wood anatomical and hydraulic traits, within‐population variation was substantial in some traits, and a boundary analysis revealed that efficient water transport was associated with vulnerable xylem and low wood density, indicating that an optimal combination of high water transport efficiency and high cavitation resistance is not possible. Our results suggest that assisted migration prescriptions may be advantageous under warming climate, but pronounced trade‐offs between survival and cold hardiness require a careful consideration of the distances of these transfers.

also experienced a reduction in precipitation in the past 25 years (Mbogga, Hamann, & Wang, 2009), and the trend toward drier conditions may continue during the 21st century (Wang, Hogg, Price, Edwards, & Williamson, 2014). As a consequence, trees may become increasingly maladapted to new climate conditions. One way of accommodating changes in climate is the use of seed sources from areas already adapted to warmer temperatures as part of regular reforestation programs. This usually implies selecting seeds from southern areas to be planted in a more northern region (Millar, Stephenson, & Stephens, 2007). Such assisted migration prescriptions depend on identifying well-adapted genotypes from matching climate regions. This can be carried out using provenance trials, in which seed sources collected from different geographic regions and different environments are planted in a common garden where genetic differences between populations may be observed. If promising genotypes can be identified in provenance trials, then these genotypes could be moved and planted where their characteristics match the anticipated climate.
Previous work on white spruce provenances in different parts of its distribution suggests that transfers toward the north can increase growth rates (e.g., Gray et al., 2016;Lesser & Parker, 2004;Li, Beaulieu, & Bousquet, 1997;Lu et al., 2014;Rweyongeza, Yang, Dhir, Barnhardt, & Hansen, 2007). A recent study proposed relatively short northward transfers for Alberta, with growth and survival of transferred seed sources putatively limited by cold temperatures in the north of the province (Gray et al., 2016). All these studies analyzed the response of tree growth to different planting environments, but there is little understanding of the physiological causes of different local adaptations.
In the boreal forest, increased temperatures may have a positive effect on tree growth, as has been observed in some white spruce populations (Danby & Hik, 2007;Lloyd & Fastie, 2002;MacDonald, Szeicz, Claricoates, & Dale, 1998). But this positive effect will only occur with adequate water availability, as the opposite effect was found in drier areas or years, showing that drought can be an important limitation for white spruce development in the future (Barber, Juday, & Finney, 2000;Chen et al., 2017;Danby & Hik, 2007;Jiang, Huang, Stadt, Comeau, & Chen, 2016;Lloyd & Fastie, 2002). Even if precipitation rates are not affected by climate change, increased temperatures will enhance drought stress in plants by increasing transpiration.
With higher transpiration, water reserves will deplete faster resulting in a heat-induced drought (Breshears et al., 2005). Moreover, snow reserves will melt earlier, further reducing water availability later in the growing season (Barnett, Adam, & Lettenmaier, 2005). Finding a productive and drought-resistant genotype might be a difficult challenge as trade-offs between growth and heat/drought resistance have been reported in white spruce (Bigras, 2000(Bigras, , 2005. The trade-off between hydraulic safety and efficiency of the xylem was analyzed in detail by Gleason et al. (2016), arriving at the conclusion that although the correlation between both traits is not always clear, the combination of both high efficiency and high resistance is not possible. This tradeoff can be partially explained by anatomical features of the tree such as wood density, conduit size, or ratio between photosynthetic and conductive tissue (Gleason et al., 2016;. Even though frost events are lower in frequency and severity under recent climate warming, extreme cold events may still occur on rare occasions, especially if overall variability in climate increases. A single unexpected frost event can cause great damage to forests if it occurs after the start of the growing season (Gu et al., 2008). Man, Kayahara, Dang, and Rice (2009) also reported severe frost damage in a white spruce stand after a late spring frost. As such frost events that cause dieback and mortality are rare, it remains difficult to assess the risk involved in moving planting stock north, even with data from long-term provenance trials because mature trees may not be as susceptible to frost damage as seedlings and saplings. Generally, differences between provenances in the onset of cold hardiness in fall are greater than in the release of cold hardiness in spring, so a movement in latitude might have a bigger effect in changing susceptibility to early frosts in fall (Aitken & Hannerz, 2001). Cold hardiness heavily relies on the phenology of the onset and release of dormancy, and a trade-off between growth and cold hardiness is usually driven by how long trees extend their growing season in the fall (Howe et al., 2003). The effect of climate change in fall usually gets less attention than other seasons even though fall events can have an important ecological impact (Gallinat, Primack, & Wagner, 2015).
While growth performance of white spruce provenances has been well studied, there is a lack of understanding of which physiological and anatomic traits are responsible for those genetic population differences. Trade-offs between growth and cold hardiness or drought resistance could pose an additional challenge for forest managers to maintain the productivity and health of our forests under climate change. The research approach of this study was to screen groups of contrasting provenances from a wide variety of climatic source environments for a broad suite of physiological and anatomical traits that are putatively adaptive. To further enhance the probability of finding trade-offs and discover genetic differentiation in adaptive traits, we selected provenances with contrasting combinations of growth and survival in a common garden field trial. The specific goals were to (1) quantify genetic population differentiation among provenances from across the entire range of the species for a wide selection of hydraulic, anatomical, and cold hardiness traits, (2) detect possible relationships among resistance to climate (cold and drought), tree growth, and survival, and (3) analyze how these traits are related to the climate of origin of the provenances. The results of this range-wide experiment could point to key traits for climate adaptation that could serve as a reference for more geographically limited studies with higher sample densities to support regional assisted migration prescriptions.

| Plant material
Plant material for this study came from two contiguous white spruce (Picea glauca [Moench] Voss) trials in a single site in central Alberta, Canada (55°17′N, 113°10′W). Both trials belong to a provenance trial experiment described in detail by Rweyongeza et al. (2007). One of the trials is planted with provenances from only the province of Alberta, while the other trial has provenances covering the whole Canadian distribution of white spruce, with some provenances present in both trials.
Both trials have five blocks each, with nine trees per block for the Alberta trial and five trees per block for the Canada-wide trial. They were established in 1982 with four-year-old seedlings. Height was assessed after 32 growing seasons in the field, and survival was calculated as the ratio of living trees to the total planted. We selected two provenances from five distinct regions that covered most of the white spruce distribution (Figure 1, Table 1). The provenance selections for Alberta corresponded to the three regions that showed genetic differences in previous studies (Rweyongeza et al., 2007): northern Alberta (nAB), central Alberta (cAB), and Foothills (FH). The other two provenances belong to climatically and geographically different parts of the species range: Yukon (YU) in the north and Ontario (ON) in the south.
The same sample trees were used for all measurements so that an individual comparison between variables was possible. A total of seven trees (at least one from each block) were selected for each provenance, except for provenances nAB.2 and cAB.2 that were present in both trials, and for which five trees from each trial and provenance were selected. The sample size was reduced for time-consuming anatomic measurements (for details, see notes in Table 2). Climate data for the provenances and study site were extracted with the software ClimateNA v5.21, available at http://tinyurl.com/ClimateNA (Wang, Hamann, Spittlehouse, & Carroll, 2016). We used the standard reference normal period of 1961-1990 as a representation of the climate of origin to which populations are putatively adapted. The period has good weather station coverage and largely precedes a strong anthropogenic warming signal.

| Hydraulic measurements
Samples for hydraulic measurements were collected between 19 May and 8 June 2015. Distal branches from sun-exposed parts of the tree crown were cut with a pole pruner. The branches were packed in plastic bags with wet paper towels to avoid desiccation and transported to a cold room (+4°C) the same day. All hydraulic measurements were made within 1 week after collection. We follow methodology described in Hacke and Jansen (2009) and Schoonmaker, Hacke, Landhausser, Lieffers, and Tyree (2010) for conductivity measurements. Briefly, we cut approximately 15-cm-long branch segments under water and attached them to a conductivity apparatus that applied a 20 mmol/L KCl+ 1 mmol/L CaCl 2 solution under controlled pressure to the segment. The outflow was measured every 10 s with an electronic balance (CP225D; Sartorius, Göttingen, Germany) until it stabilized. The average of the last five measurements was used to calculate hydraulic conductivity (K H ) with the following expression: Hydraulic conductivity was first measured to assess the native conductivity of the segments in the field. Then, samples were subjected to a partial vacuum in the measuring solution to remove any native embolism and determine the maximum conductivity. Native embolism was calculated as the ratio of the initial to the maximum conductivity.
Vulnerability to cavitation was assessed by applying a known pressure to the segment using a centrifuge and calculating the percentage loss of conductivity (PLC) relative to the maximum conductivity. This procedure was repeated for six pressure levels from −2 to −7 MPa, and the results were fitted to two commonly used functions: the Weibull function (Cai, Li, Zhang, Zhang, & Tyree, 2014) and the exponential sigmoidal function (Pammenter & Vander Willigen, 1998). The best fitted function was selected in each individual case. From these curves, we calculated the pressure at which the xylem loses 50% of its maximum conductivity (P50). Maximum K H was used to calculate xylem-areaspecific conductivity (K S ) and needle-area-specific conductivity (K L ).
Xylem area (A X ) was measured in the center of the segments using a stereomicroscope (MS5, Leica, Wetzlar, Germany). For the estimation of needle area (A L ), we first measured the projected area of a subset K H = Water flow × segment length Pressure head F I G U R E 1 Locations of origin of provenances and the common garden test site where they were grown. The area in gray delineates the species range of white spruce Ontario Foothills CentralAB NorthAB Yukon Study Site of needles, and then, the same subset was weighted after drying the needles in an oven. Using the ratio of area to weight from the subset of needles, we could estimate the total area of the needles distal to the segment by drying and weighting the remaining needles. K S and K L were then calculated according to Tyree and Zimmermann (2002): Lastly, the ratio of leaf area to xylem area (A L :A X ) was also calculated as a measure of hydraulic efficiency.

| Wood anatomy
The same segments used for hydraulic measurements were also used for wood anatomy analysis. Tracheid lumen diameter was measured for the most recent complete two rings of each segment using a radial file of three cells wide with images taken with a Leica DM3000 microscope at 20× magnification. Only noncompression wood was analyzed. The hydraulically weighted mean lumen diameter (T diam ) was calculated using the Hagen-Poiseuille formula: where d is tracheid diameter and n is the total number of tracheids measured. Tracheid length (T length ) was measured using chemical

| Cold hardiness
Healthy, sun-exposed branches were collected for cold hardiness measurements from the same trees on 22 September and brought to a cold room (+4°C) in plastic bags with wet paper towels within the same day. The whole-plant freeze testing method (Burr et al., 2001) was used to estimate frost damage. Distal parts of the branches of approximately 20 cm length were frozen at different temperatures in a programmable freezer (85-3.1A; ScienTemp, Adrian, MI, USA). We used a cooling rate of 5°C/hour and maintained the samples at the T A B L E 1 Geographic location of the provenance origins and the common garden test site used in this study Tr Length = Tracheid length (mm, N = 1 tree/provenance; 200 tracheids/tree); Tr Diam = Tracheid diameter (μm, N = 2); Density = wood density (g/cm 3 , N = 7); A L :A X = leaf-to-xylem area ratio (cm 2 /mm 2 , N = 4); Emb Nat = native embolism (%, N = 7); K S = xylem-specific maximum conductivity (mg mm −1 s −1 kPa −1 , N = 7); K L = leaf-specific maximum conductivity (mg mm −1 s −1 kPa −1 , N = 4); P50 = pressure at which 50% of the conductivity is lost (MPa, N = 7); Cold30 and Cold50 = frost damage at −30 and −50°C, respectively (%, N = 7).
test temperature for one hour, followed by warming up to room temperature at the same rate of 5°C/hr. Three test temperatures were used: −30°C, −40°C, and −50°C. After the freezing treatment, samples were transferred to a growth chamber where they were stored in transparent plastic bags with wet paper towels. Frost damage in the needles was assessed visually 2 weeks after the treatment as a percentage of discolored needle surface.

| Data analysis
All statistical analyses were conducted with the R v3.2 programming environment (R Core Team 2017). We used a mixed-effects model to calculate the means of the provenances using the variable of interest as a fixed effect and block as a random effect with the lme4 package (Bates, Mächler, Bolker, & Walker, 2015). Means were calculated as least squares means with the lsmeans package (Lenth, 2016). Multiple comparisons of the means were carried out with the general linear hypothesis method of the multcomp package (Hothorn, Bretz, & Westfall, 2008). Differences among provenances in survival were tested for statistical significance with Fischer's exact z-test and Holm's adjustment for multiple inference, implemented with the p.adjust function of the R base package. We visualized the effect of the climate of origin on the performance in the field using a multivariate indirect gradient analysis.
The ordination space consisted of standardized values of height and survival, defined as standard deviations from the overall experimental mean set to zero. We calculated the individual correlations of each climatic variable with this ordination space using the function vf() from the package ecodist (Goslee & Urban, 2007) and plotted them as vectors in the figure. To analyze trade-offs and relationships among specific hydraulic and anatomical traits, we conducted boundary analyses with quantile regressions using the function rq() from the package quantreg (Koenker, 2017). A quantile regression can be used when a response cannot change by more than some upper limit, but may change by less when other unmeasured factors are limiting (Cade & Noon, 2003). A quantile regression is similar to a linear regression, but it allows for the estimation of a quantile instead of the average. In this case, we calculated the 5% (P50) and 95% (wood density) quantiles. This regression can identify a region of a scatter plot where it is very rare to find points and can point to a trade-off even if the overall correlation between variables for the majority of the sample points is weak.

| RESULTS
Provenance field testing results conformed to the expectation that Among the physiological and wood anatomical variables measured, we found the strongest differentiation in cold hardiness ( Figure 3, Table 2). At the time of sampling, the −50°C test temperature resulted in the greatest range of observed damage from 1% to 96% in one of the Yukon (YU.2) and one of the Foothills (FH.1) provenance samples, respectively. Provenances within the same region generally behaved similarly, except for the Alberta Foothills, which showed a statistically significant difference between its provenances. Most of the provenances showed good resistance to −30°C temperatures at the time of assessment in mid-September, and only Ontario provenances showed significant damage (12%). We observed no significant differences among provenances in hydraulic and anatomic variables putatively related to adaptation to drought, except for leaf-to-xylem area ratio, due to high intrapopulation variance (Table 2).
Frost hardiness was strongly related to survival (Figure 4), but it did not show a significant trade-off with tree height (Table 3). There was a significant correlation between height and both K L and P50, but we did not find any relationship between survival and hydraulic variables (Table 3). Despite weak correlations, the boundary analysis suggested trade-offs between xylem-specific conductivity (K S ) and wood density ( Figure 5a) and between xylem-specific conductivity and the pressure at which the xylem loses 50% of its maximum conductivity (P50) (Figure 5b). Using quantile regressions, we are able to estimate a boundary (dashed line in Figure 5) beyond which values are very rare to find. Low transport efficiency (low K S values) was associated  Figure 5b). The leaf-to-xylem area ratio appeared as an influential trait for hydraulic variables as it was significantly correlated to both efficiency (K S and K L ) and cavitation resistance (P50; Table 3). Leafspecific conductivity was strongly correlated to P50, so that the more vulnerable branches had also a poorer ability to provide water to the needles. Wood density was the most influential among anatomic variables, as it was significantly correlated to tree survival and cold hardiness. Tracheid size was not related to any variable measured.
The climate of origin was strongly associated with the field performance of provenances. Climate vectors significantly correlated to height, survival, or a combination of both traits as per indirect gradient analysis. This is displayed in Figure 2, where the vector length with respect to each axis represents the strength of the association. The variables that were most strongly correlated with survival are, in this order, mean average precipitation, mean coldest month temperature, and degree-days below 0°C, suggesting that provenances with better survival came from dry climates with cold and long winters. For growth, the most important predictor variables were frost-free period and mean annual temperature at the origin of the seed source. If we evaluate the influence of the beginning and end of the frost-free period separately, we can see the higher influence of the end of growing season variable, which is also statistically significant (Table 4, Figure 6).
Annual heat-moisture index was not significantly correlated to either height or survival alone (Table 4), but it was associated with high combined survival and growth of central Alberta provenances (Figure 2).
Cold hardiness was also associated with the climate of origin.
Damage at −50°C was mostly related to day length and winter temperatures, having the highest correlations with latitude, mean coldest month temperature, and degree-days below 0°C (Figure 7a, Table 4).
Damage at −30°C was more correlated to variables corresponding to annual temperatures like evapotranspiration and mean average temperature and also to the date of the first frost event (Figure 7b,  (Figure 7b). We did not find any significant associations between climate at the origin of seed sources and hydraulic or anatomic variables (Table 4).

| Local adaptation and trade-offs
Although this study was not designed to study local optimality or to be used to derive practical assisted migration guidelines, provenances that originated near the vicinity of the study site showed the best combination of long-term performance in growth and survival in the field. Local provenances (central Alberta) originated slightly southeast of the study site which corresponds to a short northwest transfer as was recommended by Gray et al., 2016 to account for climate change that has already occurred over the last several decades. The provenances that showed both poor growth and low survival came from the Alberta Foothills (lower-left quadrant in Figure 2

| Strong population differentiation in cold hardiness
The adaptive trait with the strongest genetic population differentiation was cold hardiness, and trait values showed strong associations to latitude and the mean coldest month temperature of seed origin.
Cold hardiness was measured in branches collected around mid-September, and at this time of year, only provenances from Ontario exhibited some frost damage after exposure to -30°C. Ontario provenances usually do not experience any frost at this time of year at their location of origin. Hence, there is no need to develop cold hardiness so early in the fall when trees grow in their native climate. The other provenances do regularly experience mild frost events earlier in September in their native climate and exhibited no significant damage when exposed to -30°C at the test site. The fact that cold hardiness patterns could largely be explained by the climate of the provenance origins indicates that cold hardiness was genetically differentiated among populations, presumably mediated by day length, and that there is little potential for phenotypic plasticity in response to climate conditions. These results conform to many other studies (Howe et al., 2003 and literature cited therein). Of all the traits measured in this study, cold hardiness appeared as the most important adaptive trait, clearly linked to survival of provenances in the field.
While cold hardiness was not significantly correlated to tree height, we found a strong correlation of height with the end of the growing season at the provenance origin. This suggests that the timing of the onset of dormancy plays an important role in the growth potential of white spruce.
T A B L E 3 Pearson correlation coefficients for the mean values of provenances for anatomy, hydraulic, cold hardiness, and performance variables. Statistically significant correlations at p < .05 are highlighted in bold Cold30 and Cold50 = frost damage at −30 and −50°C, respectively, A L :A X = leaf area to xylem area, K L = leaf-specific conductivity, K S = xylem-specific conductivity, P50 = vulnerability to cavitation expressed as the pressure at which 50% of the maximum conductivity is lost, Density = wood density, Tr Diam = tracheid diameter, Tr Length = tracheid length.

| Wood anatomical and hydraulic parameters
While it seems that patterns in survival and tree height could largely be explained by cold tolerance and growing season length, respectively, we also measured a wide range of hydraulic and wood anatomical parameters, putatively related to drought resistance. Wood density emerged as a potentially influential variable as it was correlated with survival and cold resistance. Lighter wood may increase the risk of stem breakage (Spatz & Bruechert, 2000), which may be important under heavy snow loads (Hlasny, Kristek, Holusa, Trombik, & Urbancova, 2011). High wood density also provides a stronger defense against pathogens as well as lower vulnerability to drought stress (Chave et al., 2009;Hacke, Sperry, Pockman, Davis, & McCulloch, 2001;Hacke et al., 2015;Rosner et al., 2014). Higher wood densities are often found in environments that are associated with stress. Positive correlations of wood density with survival and cold hardiness seem consistent with these reports.
Tree height was negatively correlated with leaf-specific conductivity (K L ), that is, taller trees tended to show lower K L values than shorter trees. Lower values of K L corresponded to higher needle area distal to the section measured relative to its ability to transport water. K L was mostly regulated by needle area rather than by a change in hydraulic conductivity as we can see from the significant correlation between K L and A L :A X and the absence of a correlation between K L and K S (Table 3). Similar findings were reported from studies on pine populations sampled across climate gradients (Lopez, Cano, Choat, Cochard, & Gil, 2016;Martínez-Vilalta et al., 2009). At a particular transpiration rate (E), a branch with higher K L will be able to maintain a smaller water potential gradient (ΔΨ) than a branch with lower K L , because E = K L ΔΨ (Tyree & Zimmermann, 2002). While having high K L may seem conservative and advantageous from a hydraulic perspective, it comes with the disadvantage of having less photosynthetic area or having to invest more resources into nonphotosynthesizing xylem tissue to increase transport capacity.  Cold30 and Cold50 = frost damage at −30 and −50°C, respectively, A L :A X = leaf area to xylem area, K L = leaf-specific conductivity, K S = xylem-specific conductivity, P50 = vulnerability to cavitation expressed as the pressure at which 50% of the maximum conductivity is lost, Density = wood density, T Diam = tracheid diameter, T Length = tracheid length. MAT= mean annual temperature, MWMT = mean warmest month temperature, MCMT = mean coldest month temperature, TD = continentality or temperature difference between MWMT and MCMT, MAP = mean annual precipitation, MSP = May-to-September precipitation, AHM = annual heat-moisture index (MAT + 10)/(MAP/1000), SHM = summer heat-moisture index MWMT/(MSP/1000), DD<0 = degree-days below 0°C, DD>5 = degree-days above 5°C, bFFP = the day of the year on which FFP begins, eFFP = the day of the year on which FFP ends, FFP = frost-free period.
In our study, Yukon and Foothill provenances exhibited the highest K L values. Both Yukon and Foothills have the shortest frost-free periods. Trees in these regions must deal with frozen soils for the longest time, and they may also experience many freeze-thaw events in spring. Drought stress can develop in the winter when water uptake from a frozen/cold soil is impaired while needles continue to lose water (Mayr, Hacke, Schmid, Schwienbacher, & Gruber, 2006). Having a higher K L may reduce the tension in the xylem during this period and may therefore be advantageous over a greater emphasis on leaf area and growth potential.
Hydraulic conductivity standardized by xylem area (K S ) was not correlated with growth or survival, but showed complex relationships with wood density and P50. The boundary analysis approach shown in Figure 5 illustrates that high wood densities and K S values were mutually exclusive. This points to a strength versus hydraulic efficiency trade-off and is understandable given constraints arising from the geometry of tracheid-based xylem (see detailed discussion in Pittermann, Sperry, Hacke, Wheeler, & Sikkema, 2006;Hacke et al., 2015). High hydraulic conductivity (K S ) values also came at the expense of increased vulnerability to cavitation. The safety versus efficiency trade-off is complex (Gleason et al., 2016), but our boundary analysis indicates that certain trait combinations are unattainable despite high within-population variation in either of these traits.

| CONCLUSION
In this range-wide common garden study, we identified traits that could play an important role in trade-offs between growth potential and resistance to climatic stress. We should note that the growing conditions of the test site can affect what trade-offs are revealed.
Differences in adaptive traits that did not have an effect on growth and survival could be consequential under different growing environments. Wood density and leaf-specific conductivity were identified as potentially important traits, but fall cold hardiness stood out as a key trait for tree survival and also showed easily interpretable associations with climate of the provenance origin. Furthermore, the relation between height and the end of the growing season points to the importance of the timing of growth cessation in the growth potential of white spruce. Therefore, regional studies that make use of a higher sampling density to develop assisted migration guidelines should focus on measuring the effects of seed transfers on fall phenology and match provenances to new locations so that their synchronization of the onset of frost hardiness matches new environmental conditions under climate change.

ACKNOWLEDGMENTS
We are grateful to Stefan G. Schreiber for assistance with field and laboratory work. Funding was provided by an NSERC strategic grant STPGP-494071 to AH and UH, by a grant from the Faculty of Graduate Studies and Research at the University of Alberta (to UH) and by the University of Alberta Water Initiative. This study would not have been possible without generous in-kind contributions by our supporting organization, Alberta Agriculture and Forestry, with specific thanks to Deogratias Rweyongeza for logistical support and data sharing.

CONFLICT OF INTEREST
None declared.

AUTHORS' CONTRIBUTIONS
AH, UH, and JS conceived the study and designed the methodology. JS collected the data and led the writing of the manuscript. All authors analyzed the data and contributed to draft versions of the manuscript.