Structure is more important than physiology for estimating intracanopy distributions of leaf temperatures

Abstract Estimating leaf temperature distributions (LTDs) in canopies is crucial in forest ecology. Leaf temperature affects the exchange of heat, water, and gases, and it alters the performance of leaf‐dwelling species such as arthropods, including pests and invaders. LTDs provide spatial variation that may allow arthropods to thermoregulate in the face of long‐term changes in mean temperature or incidence of extreme temperatures. Yet, recording LTDs for entire canopies remains challenging. Here, we use an energy‐exchange model (RATP) to examine the relative roles of climatic, structural, and physiological factors in influencing three‐dimensional LTDs in tree canopies. A Morris sensitivity analysis of 13 parameters showed, not surprisingly, that climatic factors had the greatest overall effect on LTDs. In addition, however, structural parameters had greater effects on LTDs than did leaf physiological parameters. Our results suggest that it is possible to infer forest canopy LTDs from the LTDs measured or simulated just at the surface of the canopy cover over a reasonable range of parameter values. This conclusion suggests that remote sensing data can be used to estimate 3D patterns of temperature variation from 2D images of vegetation surface temperatures. Synthesis and applications. Estimating the effects of LTDs on natural plant–insect communities will require extending canopy models beyond their current focus on individual species or crops. These models, however, contain many parameters, and applying the models to new species or to mixed natural canopies depends on identifying the parameters that matter most. Our results suggest that canopy structural parameters are more important determinants of LTDs than are the physiological parameters that tend to receive the most empirical attention.

Canopy temperature distributions can be measured using several methods, including temperature-sensing probes such as thermocouples (Linacre, 1967), infrared thermography (Faye, Rebaudo, Yánez-Cajo, Cauvy-Fraunié, & Dangles, 2016;Leuzinger & Körner, 2007;Pincebourde & Suppo, 2016;Scherrer, Bader, & Körner, 2011), and isotopic analyses (Helliker & Richter, 2008). These studies generally have reported high variation in leaf temperatures between species (Michaletz et al., 2016 but see Helliker & Richter, 2008), and a few also illustrate within-canopy thermal heterogeneity (Leuzinger & Körner, 2007;Scherrer et al., 2011). However, all of these approaches have drawbacks that prevent them from revealing the full extent of insect-relevant variation in temperature within canopies: It is difficult to get extensive spatial coverage with thermocouples, IR thermography provides 2D views (top or side) of 3D canopies, and isotopic analyses integrate leaf conditions over large spatial and temporal scales. For estimating effects on canopy insects, we need better estimates of thermal diversity (i.e., LTDs) within individual canopies, which sets the local bounds of temperatures available to individual insect herbivores and to their local populations.
Besides mean temperature, which is of primary interest in many contexts, other aspects of variation in temperature come into play when local populations are exposed to environmental change. For example, changes in the spatial (or temporal) variance in temperature can mask or reverse the effect of a change in mean temperature (Benedetti-Cecchi, Bertocci, Vaselli, & Maggi, 2006;Dillon et al., 2016;Vasseur et al., 2014). In addition, temperature variance at local scales defines the range available for organisms (Faye et al., 2016;Pincebourde & Suppo, 2016;Woods, Dillon, & Pincebourde, 2015) to use during behavioral thermoregulation, which is a key process by which ectotherms achieve high rates of performance (Woods et al., 2015) and avoid temperature extremes (Kearney, Shine, & Porter, 2009;Sunday et al., 2014). Finally, patterns of spatial autocorrelation in temperature matter for mobile organisms .
In particular, ectotherms may pay lower energetic costs when favorable temperatures are less aggregated (i.e., are more patchily distributed) in space (Faye, Rebaudo, Carpio, Herrera, & Dangles, 2017;Sears et al., 2016). Several indices of aggregation have been developed by landscape scientists but have been used only rarely by ecologists working at local (Faye et al. 2016, Faye et al., 2017: scale of a crop field) and small scales (Caillon, Suppo, Casas, Woods, & Pincebourde, 2014: scale of a single leaf).
Both multilayer and big-leaf models, however, are effectively onedimensional models that describe changes along the vertical axis of a canopy rather than throughout full 3D canopies. Focusing on much smaller spatial scales, Saudreau et al. (2017) recently modeled how interactions between irradiance and leaf microtopography drive patterns of temperature variation within single leaves. Temperature variation at such small scales may also be exploited by individual insects (Caillon et al., 2014), but that spatial scale also does not capture the intracanopy scale that is relevant to most mobile insect herbivores.
The RATP model (Radiation Absorption, Transpiration, and Photosynthesis;Sinoquet, Le Roux, Adam, Ameglio, & Daudet, 2001) provides a good platform for bridging this spatial gap. This mechanistic biophysical model discretizes space (the unit cell is called a voxel) with a predefined dimension, allowing the user to adapt the spatial scale at which the model predicts leaf temperatures. The use of voxels makes it possible to apply the Beer-Lambert law of extinction of light across vegetation and to lower computational time considerably compared to a model simulating leaf-scale processes (see Saudreau et al., 2017). This model has been used successfully to estimate LTDs in apple trees and several other cultivated species at a voxel size of 20 cm (Ngao, Adam, & Saudreau, 2017). In addition, it has been coupled to more detailed models of heat budgets of insectbuilt structures (leaf mines) and of fruits in apple trees (Pincebourde, Sinoquet, Combes, & Casas, 2007). It has not yet been applied to trees in the wild or to forest stands, but doing so presents few conceptual problems.
The RATP model has a large number of parameters and explicitly uses environmental forcing variables that fall into three groups: (1) environmental-air temperature, wind speed, relative humidity, amounts of incoming short-and longwave irradiance, and the relative proportions of direct vs. diffuse irradiance; (2) structural-leaf area density and leaf inclination angle distribution; and (3) physiological-stomatal effects on evapotranspiration via stomatal sensitivity to temperature, water status, vapor pressure deficit, and light. It is not obvious a priori which parameters deserve the most attention.
Clearly, environmental parameters strongly influence leaf temperature (Monteith & Unsworth, 2008;Nobel, 1999), and they are among the easiest to collect. However, even the relative ordering of those variables is unclear (e.g., is LTD influenced more by variation in air temperature or incident irradiation?). Many studies have considered the impact and relative importance of canopy structure on leaf functioning (Baldocchi et al., 2002;Caldwell, Meister, Tenhunen, & Lange, 1986). However, most of these focus on photosynthesis and transpiration (Niinemets, Keenan, & Hallik, 2015). Even though leaf temperature is an integral part of many of the models used, few authors have analyzed leaf temperature explicitly, including its distribution at the scale of a single canopy or continuous canopies. This is a key gap that prevents us from connecting leaf temperature distributions to other trophic levels.
Here, we present a Morris sensitivity analysis (Campolongo, Cariboni, & Saltelli, 2007;King & Perera, 2013;Morris, 1991) of how LTDs depend on climatic, structural, and physiological parameters in the RATP model. The Morris analysis subsamples combinations of parameters from throughout the full, multidimensional parameter space and then ranks individual parameters according to their weight on the output variable (Morris, 1991). We developed our analysis using the well-described architecture and physiology of apple trees, Malus domestica (Saudreau et al., 2013). Apple is probably the plant species with the most detailed dataset on its architecture Sinoquet et al., 2009), although quite complex methodologies have been developed to capture the architecture of forest trees (e.g., Lintunen, Sievänen, Kaitaniemi, & Perttunen, 2011). We also used a common voxel size (20 cm), which was shown to predict well the leaf temperatures at local and whole canopy scales. We first apply the Morris analysis to the canopy of a single isolated tree and then to a virtual stand composed of several trees with contiguous canopies. Finally, we leverage the RATP model to provide 2D views of the 3D distribution of leaf temperatures in the canopies. These 2D views approximate how infrared devices image canopies from a distance (seeing only leaves on, e.g., the top or side). The comparison of LTDs from 2D and 3D views quantifies the error in LTDs when using remote sensing.

| RATP model
The RATP model was designed to simulate the spatial distribution of radiation and leaf-gas exchanges within a plant canopy as a function of its geometry, the surrounding climate, and the physical and physiological properties of leaves ).
The RATP model is based on many equations and requires many inputs: 10 above-canopy forcing climatic variables, more than 20 parameters related to the optical and physiological properties of the leaves, and five input parameters describing the spatial distribution of leaves within the canopy. The main model features are summarized below, and a fuller description of the model and equations can be found in Sinoquet et al. (2001) and in Appendix S1.
In the RATP model, a turbid medium approach is used to compute radiation transfer in a canopy composed of one or several species. Briefly, the model accounts for down-welling shortwave radiation from the sun (direct and diffuse) in PAR and NIR wavebands, atmospheric longwave radiation (sky irradiance), longwave radiation from the soil, and longwave exchange among leaves. Whenever possible, down-welling fluxes are set with measured values. If these are not available, they are estimated from empirical relationships based on classical meteorological data. For instance, the partition of the sun irradiance between direct and diffuse parts can be inferred from the clear sky index (Reindl, Beckman, & Duffie, 1990), and the sky longwave irradiance from air temperature and vapor pressure of water in the air (Iziomon, Mayer, & Matzarakis, 2003). The plant canopy is embedded in an array of 3D voxels (set to the optimal dimensions of 20 × 20 × 20 cm; Sinoquet, Sonohat, Phattaralerphong, & Godin, 2005), and the canopy structure is entirely defined by the spatial distribution of the leaf surface area and the inclination angles of the leaves. Each individual voxel is characterized by a leaf area density and composed of a sunlit leaf surface area that intercepts both direct and diffuse radiation and a shaded leaf surface area that intercepts diffuse radiation only.
For simplicity, we considered canopies composed of leaves only, and all other organs such as branches and fruits were neglected.
To estimate local leaf temperature, the energy balance equation between incoming and outgoing fluxes is closed in each voxel for sunlit and shaded leaves.
Achieving energy balance requires additional information about sensible and latent heat fluxes (convective transfer between leaf and air, and transpiration cooling, respectively). In the RATP model, sensible heat flux lost or gained by a leaf is assumed to be related to a boundary layer conductance and the difference in temperature between the leaf and the surrounding air (Monteith & Unsworth, 2008).
The boundary layer conductance is linearly related to the local wind velocity (Daudet, Le Roux, Sinoquet, & Adam, 1999). We accounted for wind attenuation in the canopy using an empirical relationship between relative wind speed and the cumulative leaf area along the wind path (Daudet et al., 1999). The transpiration rate of a leaf is driven by the physiological response of the plant (stomatal conductance) and the evaporative demand (Monteith & Unsworth, 2008). Stomatal conductance is controlled by the microclimate at the leaf surface (PAR intercepted, CO 2 and VPD), the leaf temperature, and the leaf physiological state following the Jarvis model (Jarvis, 1976) assuming no interactions between these variables on the stomatal response. The overall leaf conductance for water vapor transfer is computed by combining boundary layer and stomatal conductance of both upper and lower leaf surfaces. The overall set of equations in the RATP model consists of a nonlinear system with unknown variables, the temperature of sunlit or shaded leaf areas, to be solved (Appendix S1). Brent's iterative method (Brent, 1974) is used to solve the model. Model outputs of concern are the irradiance at the leaf surface and temperature of the shaded and sunlit area of foliage at the voxel scale.
In all subsequent analyses, we analyzed patterns of thermal variation in the entire 3D canopy of the forest stand (results for single tree are presented in Supporting Information) and just in the 2D views from above. By comparing 2D and 3D patterns of variation, we were able to estimate the utility of flyover views (e.g., what an infrared-equipped UAV would see) for estimating spatial patterns of thermal variation throughout the entire 3D canopy.
F I G U R E 1 Three-dimensional view of the isolated apple tree (a) and the virtual stand (b) used in the simulations. The virtual stand was made artificially by putting together twenty isolated tree canopies to construct a continuous canopy. The total leaf surface area was 24.3 m 2 (with 10,214 leaves) and 202 m 2 (with 93,992 leaves) for the isolated tree and the continuous canopy, respectively. Space was discretized into 574 and 5,340 voxels (20 × 20 × 20 cm) for the isolated tree and the continuous canopy, respectively. A total of 107 and 623 voxels can be seen from the top of the isolated tree and the continuous canopy, respectively

| Parameter values
Because the model has been applied most extensively to apple trees, we used as the base set of parameters those that have been meas- In the stand, trees were regularly spaced on a grid of 125 cm by 150 cm. These architectural traits directly defined the LAD and LIAD parameters (Table 1). In addition, physiological parameters related to stomatal responses (Jarvis functions) were measured on these same trees (Table 1; Pincebourde et al., 2007;Saudreau et al., 2013). For the parameters in the Morris sensitivity analysis (Table 1), we took ranges of values that are typical for temperate regions where apple trees are cultivated in France.
During the simulations, optical and physiological parameters other than those used in the Morris analysis (Table 1), including the spatial arrangement of 3D voxels, were kept constant. The 3D space was discretized into voxels of 20 × 20 × 20 cm size, leading to 574 and 5,340 voxels for the isolated tree and the forest stand, respectively.
The simulations were performed for a given day (day of year 234), time (12 h GMT), and location on the earth's surface (France: 45° north latitude, and 0° east longitude). The sun was thus positioned in the sky at an elevation angle of 68° and an azimuthal angle from north of −178°.

| Statistical summaries
The statistical analysis focused primarily on the sunlit portion of the foliage, which contained most of the temperature variation, but included some analysis of shaded portions too. Shaded portions were always close to ambient air temperature. For each simulation run, several key statistical parameters were computed.
Heterogeneity can be described in terms of composition and configuration. Composition was described by the mean and variance of temperature at the canopy scale. Configuration was quantified using two spatial correlation indices: Moran's I (I Moran ) and Geary's C (C Geary ) indexes; results for Moran's I are presented in the main text and for Geary's C in Appendices S2 and S3. These indices reflect the extent to which patches (voxels) with similar temperatures are aggregated in space (Jumars, Thistle, & Jones, 1977), here at the scale of the canopy. Moran's I is based on a crossed-product centered to the overall mean, while Geary's C is sensitive to the deviation between pairs of points independent of the mean (Fortin, Drapeau, & Legendre, 1989;Jumars et al., 1977). I Moran is zero when the spatial distribution is random, and it moves toward −1 and +1 for dispersed and clustered distribution, respectively. C Geary departs from 1 (no spatial autocorrelation) down to zero or above 1 for increasingly negative or positive spatial autocorrelation, respectively.
Briefly, Moran's I detects aggregation due to extreme temperature values in several adjacent voxels, whereas Geary's C tests whether adjacent voxels contain similar temperatures (Jumars et al., 1977).
Moran's I is often seen as a global, and Geary's C as a local, index of aggregation.

| Sensitivity analyses
The Morris approach (Morris, 1991) provides an efficient way of exploring large parameter spaces. In our case, we used 13 parameters and variables with four levels each, or 4 13 = 67,108,864 unique combinations of values. Because simulating each set takes several minutes on our computers, exploring the total factorial space is effectively impossible. The Morris approach cuts down on total coverage by taking a kind of random but well-distributed walk through the 13-dimensional hypercube of parameter values. In the analyses presented below, for example, we explored only 5,600 unique parameter combinations-which is <0.01% of the total number-but were still able to identify the relative influence of individual parameters on simulation outcomes. For each of the i parameters (among environmental, structural, and physiological groups), the analysis calculates two parameters, μ * i and σ i , which refer, respectively, to the overall strength of the effect of the ith parameter on simulation output values and the strength of its interactions with other parameters (King & Perera, 2013). This analysis was applied to the different statistical outputs explained above: mean leaf temperature, canopy scale variance in leaf temperature, and the Moran and Geary indices.
Once the relative importance of parameters was established using the Morris approach, we used standard sensitivity analyses to visualize

VPD max
Stomatal response to  Table 1).

| Ranking the relative importance of parameters
The Morris analysis ranked the 13 parameters according to their μ* (overall relative impact of each parameter) and σ (the extent to which a given parameter interacts with others). Overall, the analysis indicated that the physiological parameters affected LTDs much less than did structural parameters and that the climatic parameters were the most influential ( Figure 3).
Mean leaf temperature was driven mainly by climatic parameters, and in particular air temperature (highest μ*) and wind speed (highest σ) (Figure 3). These two parameters affected mean leaf temperature similarly for top-viewed and all voxels (Figure 3a,b). The variance in leaf temperature was driven mostly by wind speed and to a lesser extent by the other climatic parameters, except humidity ( Figure 3). However, the structural parameter LAD was also important, especially for the thermal variance in the top-viewed voxel (Figure 3c,d). In general, we observed a positive relationship between the impact (μ*) and the interaction strength (σ) of parameters, except for air temperature, which had a high impact on mean leaf temperature but interacted little with other parameters.
The Morris analysis showed further that the spatial autocorrelation of leaf temperatures (I Moran and C Geary ) was affected primarily by the structural parameter LAD and by climatic variables such as air temperature, irradiance, and Kt (Figure 4). The low σ of LAD for the Moran index indicates that this structural parameter does not interact strongly with the other parameters. Again, the physiological parameters were not important for the spatial configuration of LTDs for all and top voxels (Figure 4).

| Standard sensitivity analysis
The LTDs for top and all voxels were compared by visually inspecting the quantitative difference between the respective box plots A similar trend was found for the variance in the LTDs: low wind speed and high irradiance generated large deviations between top and all voxels, and all voxels had a higher thermal variance than top voxels under these two specific circumstances (Figure 5b,d). Both LAD and Gsmax did not cause particular deviations between top and all voxels (Figure 5f,h).

| Comparing the LTDs of the canopy surface and the whole canopy
Using the mean and variance of top voxels as predictors, we were able to estimate the mean and variance of all voxels in the continuous canopy with fairly high confidence ( Figure 6). Not surprisingly, top voxels overestimated mean leaf temperatures and underestimated thermal variance, throughout the entire canopy ( Figure 6).

| D ISCUSS I ON
Plant canopies are complex, heterogeneous spaces (Finnigan, 2000;Sinoquet et al., 2001;Urban et al., 2012). Leaf temperature distributions (LTDs) in plant canopies strongly influence local distributions and performance of small organisms, including phytopathogens, insect herbivores, and their predators and parasitoids (Bernard et al., 2013;Chelle, 2005;Pincebourde & Woods, 2012;Pincebourde et al., 2007), and LTDs can play a role in the responses of these organisms to environmental changes (Pincebourde, Murdock, Vickers, & Sears, 2016;Pincebourde & Suppo, 2016). Our sensitivity analysis of the RATP model indicates, not surprisingly, that climatic parameters have strong effects on mean leaf temperature at the canopy scale. This is because leaves are strongly coupled to atmospheric conditions, in particular to air temperature. However, we also found that structural parameters, especially leaf area density (LAD), can influence both the mean and spatial variance of temperatures, in some cases being more important than climatic parameters such as relative humidity. Climate parameters also influenced variance in temperature.
Intracanopy variance in temperature is directly generated by irradiance interacting with structural aspects of the canopy (but see below). At low irradiance, variance falls to low levels; at high irradiance, variance is magnified by irradiance attenuation through the canopy (Ngao et al., 2017). The ratio of diffuse to direct irradiance (parameter Kt) also affects the variance of canopy LTDs, but this effect is much less studied compared to others (Urban et al., 2012).
More diffuse irradiance homogenizes leaf temperatures because shaded leaves receive relatively more radiative energy. Moreover, an increase in the diffuse portion is normally associated with cloud cover, which implies that the leaves exposed to the (overcast) sun are receiving less radiation at the same time (i.e., there is convergence of incoming irradiance for shaded and sunlit leaves) (Urban et al., 2012). Other climatic parameters may also influence the thermal variance. Wind speed generally is higher at the surface of canopies than inside, and this process is included in the RATP biophysical model , although turbulent regimes in forest canopies can generate complex flow (Finnigan, 2000). The wind effect depends on the density of foliage, which explains why σ is high for this factor compared to others when looking at the thermal variance. Importantly, structural parameters (especially leaf area density, LAD) were, in most cases, more important than physiological parameters. Their effects were strongest on the variance and autocorrelation indices of LTDs, meaning that architectural traits strongly influence the magnitude and configuration of spatial heterogeneity in leaf temperatures across canopy forest covers (see also Ref. (Scherrer et al., 2011)). Here too, this result emerges from interactions between irradiance and structure. Greater LAD means that much of the incoming irradiance is intercepted at the top of the tree, generating larger gradients and more aggregated temperatures. In addition, the position of the sun in the sky relative to the orientation and inclination of the leaf surfaces determines how much radiative energy is intercepted by these surfaces (geometric relationshipssee Ref. (Oke, 2002)). Therefore, a high diversity of leaf angles should produce a high diversity of leaf temperatures spread throughout the canopy. Surprisingly, however, the effects of structural parameters exceeded the influence of physiological parameters. Although transpiration rate is well known to depress leaf temperatures (Campbell & Norman, 1998;Potter, Davidowitz, & Woods, 2009)  , an approach that necessarily homogenizes variables within voxels and thereby ignores variation at finer scales (Caillon et al., 2014;Saudreau et al., 2017). Discretization is necessary however to lower computational time and to allow the application of Beer's law to calculate the attenuation coefficient of irradiance across the forest canopy . In addition, our aim was to analyze LTD at the scale of canopies, not at the scale of individual leaves (for leaf-scale variation, see Saudreau et al., 2017), although the two scales are inter-related and can be combined (Pincebourde & Suppo, 2016). Remote sensing and most plant-atmosphere exchange models are based on local/canopy scale processes, and our aim was to quantify the relative influence of the main climatic, structural, and physiological parameters in the same context as do these approaches.
Second, the Morris approach is designed to estimate main effects of (μ*) and interactions among (σ) all parameters across the entire parameter space, and it does not account for potential correlations among variables that may occur in the real world (e.g., high humidity often co-occurs with low irradiance). Thus, a more focused sensitivity analysis, in a subset of parameter space using predetermined correlation structures, may provide a different ranking of the more influential parameters. The advantage of the Morris analysis is that it provides a global view by examining the entire parameter space.
The main consequence of the interplay among canopy parameters, and of their relative influence on canopy LTDs, is that the temperatures measured at the surface of canopies are a reasonable proxy of the LTD within the entire continuous canopy. Our results indicate that over a large range of parameter space, the LTDs of topviewed voxels approximate the LTDs of all voxels. Nevertheless, low wind speeds and high irradiance both result in strong deviations between the LTDs of top-viewed voxels and all voxels, because under these conditions, the temperature of leaves at the top of the tree canopies increases markedly relative to other leaves. Still, we estimate that the prediction error of the median leaf temperature is <3°C under these extreme conditions (but see Figure 6). These deviations are reflected both in mean leaf temperature and, to a lesser extent, the thermal variance. Overall, the least influential parameters from the Morris analysis (e.g., stomatal conductance) did not induce deviation between top voxels and all voxels.
Our results support several recommendations for biologists interested in thermal ecology and the global change biology of forest insects. First, we need more data on plant architectures (Barthélémy & Caraglio, 2007) across latitudes, altitudes, and biomes. Currently, detailed architectural traits are available primarily from crops and orchards Sinoquet et al., 2009) but not from single species or mixed-species stands in the wild (see Ref. (Farque, Sinoquet, & Colin, 2001) for an exception).
Currently, therefore, it is not possible to compare our results on apple with similar approaches on forest trees. The well-known TRY database (Kattge et al., 2011), which lists numerous plant traits for a large number of species across the world, contains architectural traits, but they are underrepresented compared to physiological parameters. Several tools exist to measure basic architectural traits, including electromagnetic digitizers (Godin, Costes, & Sinoquet, 1999;Sinoquet, Thanisawanyangkura, Mabrouk, & Kasemsap, 1998), LIDAR scanning (Béland, Widlowski, & Fournier, 2014), and orthophotography for simple canopies (Rich, 1990). We therefore call for additional use of these tools to collect architectural traits on diverse groups of plant species. Second, the overall range of surface temperatures (3D) available to small organisms living on leaf surfaces can be approximated from infrared thermography of canopy surfaces (2D). Remote sensing satellites, unmanned aerial vehicles (UAVs), or crane towers equipped with infrared sensors are now providing such 2D temperature data for a wide range of plant species across latitudes (Dong, Prentice, Harrison, Song, & Zhang, 2017;Faye et al., 2017;Leuzinger & Körner, 2007). Our results suggest that remote sensing data at adequate spatial resolutions can be used to predict the distribution and performance of plant-associated organisms within forest canopies.

ACK N OWLED G M ENTS
We thank INRA, Université Clermont Auvergne, the Institut de Recherche sur la Biologie de l'Insecte, PICS (CNRS project MEGALEAF to SP, HAW, and MS), MPG Ranch (Montana), and The University of Montana for supporting this work.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
HAW, MS, and SP conceived of the project; MS parameterized and ran the RATP model and Morris analyses; HAW did additional analysis; HAW and MS constructed the figures; and HAW, MS, and SP wrote the manuscript.

DATA ACCE SS I B I LIT Y
Data from the article will be made available in an online repository once the article is accepted.