Comparison of water-use characteristics of tropical tree saplings with implications for forest restoration

Tropical forests are experiencing reduced productivity and will need restoration with suitable species. Knowledge of species-specific responses to changing environments during early stage can help identify the appropriate species for sustainable planting. Hence, we investigated the variability in whole-tree canopy conductance and transpiration (Gt and EL) in potted saplings of common urban species in Thailand, viz., Pterocarpus indicus, Lagerstroemia speciosa, and Swietenia macrophylla, across wet and dry seasons in 2017–2018. Using a Bayesian modeling framework, Gt and EL were estimated from sap flux density, informed by the soil, atmospheric and tree measurements. Subsequently, we evaluated their variations with changing vapor pressure deficit (VPD) and soil moisture across timescales and season. We found that Gt and EL were higher and highly variable in L. speciosa across seasons than S. macrophylla and P. indicus. Our results implied that water-use in these species was sensitive to seasonal VPD. L. speciosa may be suitable under future climate variability, given its higher Gt and EL across atmospheric and soil moisture conditions. With their lower Gt and EL, P. indicus and S. macrophylla may photosynthesize throughout the year, maintaining their stomatal opening even under high VPD. These findings benefit reforestation and reclamation programs of degraded lands.

Tropical forests are one of the largest carbon sinks in terrestrial ecosystems 1 . Climate change has been predicted to result in more intense with protracted dry or wet spells, with varying duration, intensity, frequency, and spatiotemporal spread in different parts of the world 2,3 . In fact, drought-induced mortality of global tropical forests has been recently documented 4 . Species mortality caused by anthropogenic factors can be important 5 and can result in a loss of carbon storage, leading to a cyclical downturn of increasing temperatures and further mortality. Therefore, many countries are taking preemptive action to reduce such impacts through reforestation. Reforestation includes planting trees on previously degraded lands, which can be demanding on newly planted trees. Under the predicted intensification of climate change, reforested trees may experience harsh conditions, resulting in an unsuccessful establishment. To ensure successful establishment in degraded lands, species robustness to seasonal variations, especially under a climate change scenario is needed. Because trees in urban conditions likely experience extreme environmental stresses from heat, drought, and flood, another interesting approach is to study the tolerance and resilience of urban tree species to climate variability which will infer the physiological responses of trees and their vulnerability to climate extremes 5,6 .
In Bangkok, street tree species differ in leaf phenology (i.e., deciduous or evergreen), habitat (ranging from hill evergreen forest to mangrove forest), and origin (native or exotic) 7 , with the most common street tree species being Pterocarpus indicus, Lagerstroemia speciosa, and Swietenia macrophylla 8 . P. indicus is a fast-growing medium-size tree native to Asia and can be either deciduous or facultative evergreen, depending on the availability of soil moisture and the openness of the growth area. L. speciosa is a fast-growing medium-size tree native to Southeast Asia and is a deciduous species, usually found in mixed deciduous forests. S. macrophylla is an

Results
Model performance. An initial sensitivity analysis was done to determine the capacitance parameter, β (i.e., β was varied between 0.22 and 1 or a storage time between 120 to 0 min), as indicated by high R 2 and significant p values. The β parameter represents the water storage capacity in stems which can be discharged through transpiration, buffering its daily fluctuations. In the context of our analysis, the response in sap flux density can be dampened by lags in stomatal response and capacitance (i.e., water storage). Additionally, we also ensured that model parameters converged to clear posterior distributions. These were obtained at different values of β for the three species, during a given season. Additionally, based on our analysis, keeping the prior on β weak (i.e., non-informative prior) improved the model predictions further. Hence, for a given species and season, we initialized the simulation with a given β value, which stabilized to a final value within the first 1000 iterations. Table 1 indicates the values with which the simulation was initialized (as determined through the sensitivity analysis) and the mean value obtained after removing the burn-in. For P. indicus and L. speciosa during the wet season, the variations in sap flux density were best explained by storage time which was around half of that during the dry season, while S. macrophylla had a higher storage time in the wet season. A lower storage time means   Table 2, whole-tree canopy conductance varied the greatest in L. speciosa with the relative sensitivity to natural log of VPD being low (λ/G ref < 0.50) during both seasons for all the species, with only exception being P. indicus during the dry season (λ/G ref > 0.70). Both G t and E L showed similar seasonal variation for a given species but with varying intensity (Fig. 2). The uncertainty in daily G t and E L averages was high during the wet season, indicating that the variability in VPD and soil moisture contributed greatly, relative to the model error. Such errors can be amplified further at much shorter, half-hourly timescale, as indicated by the variable R 2 for the three species (Fig. 1). Overall, the mean values of both G t and E L were statistically different among species and between seasons. L. speciosa had the highest G t and E L while P. indicus had the lowest values. Both G t and E L were significantly higher in wet season than in dry season. This difference was more pronounced in G t , where, for L. speciosa, it was about 10 times higher than P. indicus, while E L was approximately 2 times higher. However, only the mean daily values indicated a significant interaction between species and seasons. This may be attributed to the lower Table 2. Mean stomatal sensitivity along with mean daily, weekly and monthly G t and E L (± SD) for the three investigated species during the wet and the dry season. The same alphabetical superscript indicates no significant difference from each other and vice versa at 5% significant level. P-values in bold and having ** indicate a statistically significant difference between means of the main factors and the interaction between the main factors.  www.nature.com/scientificreports/ number of samples in the weekly and monthly analysis, leading to a large variation and undetected difference in the ANOVA test. Daily, weekly, and monthly G t and E L were the greatest in L. speciosa in either of the season and were the lowest for P. indicus ( Table 2). The mean G t and E L for L. speciosa during the wet season was around 77 mmol m −2 s −1 and 0.76 mmol m −2 s −1 , respectively and reached their maximum value (around 140 mmol m -2 s -1 and 1.6 mmol m −2 s −1 , respectively) around mid-July, 2017. This can be attributed to the low total leaf area of L. speciosa (averaged 2.55 m 2 ) compared to P. indicus and S. macrophylla, which had similar leaf areas (averaged ~ 3 m 2 ). In general, G t declined drastically for L. speciosa, with the other two species experiencing marginal declines as time progressed from wet to dry season, till the end of the study period.
Diurnal variations. During both seasons, VPD variations over the day were between 0.7 and 2.0 kPa, with VPD being significantly higher during the wet season than in the dry season from 8:00 to 17:00 local time (LT) (Fig. 3a). The soil moisture content measured in each species was lower during the dry season (Fig. 3b). The course of hourly G t and E L (averages from 6:00 to 18:00 LT) are plotted for the three species ( Fig. 3c and d, respectively). A uni-modal pattern was observed in both G t and E L in the three species, showing peaks between 8:00 and 14:00 LT ( Fig. 3c and d). Maximum G t was reached early in the day (around 9:00 LT) for L. speciosa, while the E L rates peaked approximately 2 h later during the wet season. A sudden rise to high values during the early hours was more pronounced during the wet season. After this, both parameters declined relatively consistently, during both seasons. Overall, G t and E L were consistently higher during the wet season, with the estimated E L values for L. speciosa being twice that of the dry season. A plateau during the wet season (less pronounced in the dry season), was observed for S. macrophylla (from 9:00 to 13:00 LT).
Over various time scales (weekly and monthly not shown), E L was positively correlated with G t within both wet and dry seasons ( Fig. 4 showing daily estimates only). The significant relationships revealed that, on average over both seasons, a unit change in G t caused a change of 0.45, 1.16, and 0.76 units in E L for P. indicus, L. speciosa, and S. macrophylla, respectively.
Atmospheric and soil demand. Responses of daily E L to atmospheric evaporative demand, represented by VPD, and soil moisture, were different between the two seasons and species (Fig. 5). Generally, E L increased with increasing VPD (Fig. 5a-c) in all species. The sensitivity of E L to VPD was relatively higher during the wet season and highly pronounced for L. speciosa and S. macrophylla (compare green and red lines in Fig. 5b and c), while that for P. indicus was similar in both seasons (Fig. 5a). Strong decreases in E L with increasing soil moisture was observed during the wet season for L. speciosa and S. macrophylla, while the effect was highly suppressed in P. indicus (green lines in Fig. 5d-f). The observed trends of decreasing E L with increasing soil water indicated possible flooding condition during the wet season. Some extra rainfalls may contribute to this, although we had no access to rainfall data at the site to confirm it. However, the range of soil moisture during the wet season was relatively small with differences between the maximum and the minimum around 0.1 m 3 m −3 for both species (Fig. 5e, f) and the significant trends may be driven by the 'extreme' value such as the data point at the lowest soil moisture of S. macrophylla (the leftmost green symbol in Fig. 5f). During the dry season, E L did not change significantly for both L. speciosa and S. macrophylla but it increased marginally with increasing soil moisture for P. indicus.

Discussion
Model performance to predict sap flux density in saplings. Using the Bayesian modeling technique 31 , we estimated G t and E L from sap flux density measured in saplings of three common urban species in Bangkok. The sap flux rates were closely related to the weather variations, particularly VPD. Previous studies have indicated that the stomatal response is sensitive to VPD as it directly affects the leaf water loss which, together with plant hydraulic system, controls leaf water potential and turgor pressure [39][40][41][42] . The best model performance was obtained by uniquely parameterizing each species and season with variable dampening effects (as indicated by the β parameters in Table 1), based on the dependence of water storage and hence transpiration on changes in soil and plant water status 43 . Given many uncertainties related to measurement errors (e.g., those in measuring environmental data), the overall agreement between the model estimates and the data was sufficiently accurate, as indicated by the respective R 2 and RMSE values, implying that the model estimates of G t and E L were reliable. Furthermore, VPD and soil moisture could explain a large portion of the variability in the sap flux density of the three species, as indicated by mean R 2 during wet (dry) seasons; 0.64 (0.65) for P. indicus, 0.70 (0.57) for L. speciosa, and 0.81 (0.78).
While the model output gives reasonable values of sap flux density which allows good estimates of G t and E L across the three species, certain important caveats need to be considered. In our study, assumptions were made while predicting J t , G t , and E L from the sensor data, soil moisture, and environmental variables. First, we assumed that there were no radial or azimuthal variations in sap flux density, due to the small stem size. Additionally, as the non-conducting wood (pith) was negligible, it was assumed that the sapwood covered the entire stem. Thus, the estimated values are for potted saplings and more mature individuals could have a relatively modest variation between the estimated values of G t and E L , due to a higher stomatal sensitivity to VPD 44,45 . The difference between present study and the previous study 31 that employed the StaCC model is that their estimations were based on sap flux density in mature trees and were closely related to the total leaf area index of the crown, without the need to individually quantify the leaves that constitute it. Also, due to many unusable data points, we did not use the PAR to generate the estimates. However, we tested the model run with the data provided by the original model 31   www.nature.com/scientificreports/ may thus hide errors associated with such assumptions on canopy structure and missing information through appropriate parameter values. As indicated by the R 2 and RMSE values, the model G t and E L estimates could be generated without a substantial loss in model performance.
Canopy conductance and transpiration. The distribution of G t was observed to be much broader for L.
speciosa, relative to P. indicus and S. macrophylla, and approaching maximum G t during the wet season (Table 2  www.nature.com/scientificreports/ and Fig. 2). Species-specific stomatal response is sensitive to environmental drivers, including atmospheric and soil moisture conditions, in addition to species phenology, which may promote a greater water flux during specific periods. The observation of relatively high levels of G t in L. speciosa (deciduous species 7 ) may be attributed to the fact that drought deciduous species tend to maximize carbon gain during the growing season (wet season), at the cost of a shorter leaf longevity compared to dry evergreen species 46 . P. indicus and S. macrophylla, which are evergreen species, had relatively lower values of G t and E L when compared to the deciduous L. speciosa, during the study period. The evergreen species are able to photosynthesize all year round 47 , without losing all the leaves in the canopy during any season 48 . Therefore, opened stomata for photosynthesis can be maintained at the same rate throughout the seasons (Table 2). During the wet period, we found that the atmospheric demand (VPD) was higher than in the dry period from 10 am to 6 pm (Fig. 3a) resulting in a higher canopy conductance and transpiration. Typically, VPD during the dry period should be higher than the wet period. However, in our experiments, saplings were placed on the roof top of a building, which could have had some confounding effect of heat dissipation from the building or the cement floor resulting in a high VPD during the wet season.
In the present study, as the saplings were potted, isolated, well ventilated, and placed on the balcony of a building, the variability in G t was closely related to the atmospheric demand, with an approximate but weak linear relationship observed between daily variability in G t and E L (Fig. 5) 49 . During the dry season, G t and E L were greatly reduced in all three species, which could indicate a moderation in the internal water potential, with the most prominent decrease in L. speciosa. The diurnal course of G t and E L was also variable between the three species (Fig. 5). In a previous study 50 , it was reported that the sap flux density peaked just before noon and then decreased gradually throughout the late afternoon. This could be explained by a relatively higher stomatal closure in all species under high evaporative demand.
Soil and atmospheric demand. Stomatal response to VPD is usually optimized by minimizing water loss for a given amount of carbon assimilation 51 . Apart from sustaining a high rate of assimilation, this adaptability also helps in maintaining optimum leaf temperatures 52 . A reduction in the available soil moisture can lead to stomatal closure, reducing the evaporative demand on the internal water potential and decreasing the rate of depletion in root zone water 53 . Similar relationships with soil moisture have been observed in the Amazonian 54 and Bornean rain forests 55 and may reflect a response mechanism to counteract short-term water stress induced by partial drying of root zone during dry periods.
A previous study 56 reported that S. macrophylla was able to extract water at relatively lower levels of soil water content. However, our study was on potted saplings, whose root expansion was limited in a container, where S. macrophylla was able to extract most water in the container until the soil moisture reached 0.2 m 3 m −3 in dry period (Fig. 5b). The highest transpiration rates estimated for L. speciosa might indicate a low xylem resistance to water flow and a deep root system 57 to maintain the water flow during dry season. Soil moisture used by L. speciosa during the dry season was also the highest among the three species and was above 0.25 m 3 m −3 even during times when the estimated transpiration rate was the highest. This could be due to the total leaf area of L. speciosa being the lowest, such that the rate of water depletion in the container was lower than that for S. macrophylla and P. indicus.
Suitable species for restoration under a future climate change scenario. Tropical species may be susceptible to variation in temperature, with a relatively small change in mean temperature causing a disproportionately large change in plant function due to heat related stresses. Additionally, extreme weather events due to erratic precipitation patterns can lead to an increased likelihood of excessive precipitation or protracted dry spells. The adaptability of a species to respond to such occurrences would strongly impact their survival. This trade-off can cause susceptibility to drought conditions in those plants that have been growing in waterlogged conditions and vice versa 58 . Such large variations in temporal difference in soil moisture conditions may be detrimental for endangered species having narrow physiological tolerances 59 or poor competitors 60 .
Forests degraded due to prolonged mono-crop practice, from fires, or a high mortality rate resulting from climate change, can be restored back from a state of being totally razed. Buffer zones around national parks fall under a sensitive land use criterion, and are treated as a protected area. But native people residing in such areas are allowed to live of the land while being under strict regulatory supervision. Illegal deforestation and ill-effects of climate change have resulted in a high mortality rate of forest species. One strategy for devising a restoration plan is to choose native species because the introduced alien species may easily spread and are harmful to the native species.
Deciduous species exhibit higher levels of transpiration and greater cooling and drought avoidance through leaf loss, relative to evergreen species, which have a limited physiological and morphological adaptability to survive wider swings in water and temperature ranges 61,62 . While evergreen species, mostly native to wetter environments, can have fine roots concentrated near the soil surface during the dry season as nutrients and water are readily available on the soil surface. Therefore, when water is not a limiting factor on the forest floor, evergreen species compete for nitrogen and light resources 63 , which under a restoration scenario, apart from soil water and nutrients are key limiting factors determining plant survival.
All three species in this study are native to tropical regions and can be considered for establishment under a restoration program. Moreover, since all of them are fast-growing species, they can also be pioneer species in degraded or barren lands as they would provide shade for other species. Under circumstances of low soil nutrition and the need for high rate of nutrient cycling, L. speciosa could be a suitable choice for restoration as it tends to photosynthesize rapidly, as indicated by higher canopy conductance and transpiration ( www.nature.com/scientificreports/ Fig. 2). Furthermore, given its deciduous nature, it would shed leaves during the dry season resulting in higher return rates of soil nutrient recycling. Under a climate change scenario with higher temperatures and drier atmosphere and soil water stress, P. indicus is likely to be suitable choice under the restoration program, as it maintained a steady canopy conductance and transpiration during both wet and dry season soil moisture and air vapor demand (Fig. 3). This suggests that P. indicus would maintain a steady growth rate, independent of how dry the soil or air is. In other words, the stomata of P. indicus are relatively less sensitive to air and soil moisture demands when compared to L. speciosa and S. macrophylla. Also, P. indicus would be able to survive even if the future climate is predicted to be wetter than at present. However, the only downside of considering P. indicus is its potentially relatively slow growth rate due to low stomatal opening for photosynthesis but if the annual growth is to be accounted for, it may not be much different from the other two species.

Conclusions
We report the model estimates of G t and E L for the three potted sapling tree species across temporal scales during the wet and dry season using the StaCC model, which has been previously used only with the mature tree species. With our unique parameterization for species and season, we concluded that the model was able to perform well with the sapling data and could be used further to provide meaningful estimates of G t and E L for saplings. As the estimated transpiration rates of the studied species were variable, it is possible to select the most suitable species to be planted accordingly to restore degraded lands for effective reforestation. Transpiration per unit leaf area was the highest in L. speciosa while that in S. macrophylla and P. indicus was lower in both seasons. During the wet season with high VPD, higher amount of water loss through transpiration occurred compared to the dry season with a lower VPD. P. indicus is likely to be insensitive to changes of soil moisture and atmospheric demands, while L. speciosa responded to even small variations in the weather conditions. Leaf phenology might play a crucial role in G t and E L due to differences in the species' behavior. Being a deciduous species, L. speciosa would have a stronger stomatal control, while P. indicus and S. macrophylla, as evergreen species, tended to have low water use, given the ability to photosynthesize all year round, and would be able to maintain their stomatal opening even under the high atmospheric demand. Nevertheless, future studies on field-grown trees should be performed to confirm the findings.

Methods
Experimental setup and environmental measurement. The study site was established on the balcony of the 4th floor of a building in Chulalongkorn University, Bangkok (13° 44′ 2.9′′ N 100° 31′ 54.1′′ E). Based on the 30-year-record (1981-2010) of climate data from a Bangkok metropolis station (Thai Meteorological Department), mean annual air temperature was 28.6 °C and mean annual precipitation was 1648 mm. This study was undertaken from 25th July 2017 to 11th February 2018, which included both the wet (25th July-31st October 2017) and dry season (1st November 2017 to 11th February 2018). Ten trees for each species were bought from a nursery (Chatuchak market) and an automatic irrigation system was used to water the trees once a day. Despite the daily irrigation, the fluctuations of soil moisture may be large enough to affect any physiological responses of the potted saplings, hence we considered soil moisture as a covariate in this study. A detailed description of the experimental setup and the measurements can be found in a previous study 50 which mainly reported the observed sap flux density without scaling up to canopy conductance. We measured half-hourly air temperature (°C) and relative humidity (RH in %) with a probe (HC23-L, Campbell Scientific, Logan UT, USA), which were used to calculate VPD (kPa) 64 . Photosynthetically active radiation (PAR, umol m −2 s −1 ) was measured using a quantum sensor (LI190R-L, Li-COR Biosciences, Lincoln, NE, USA. All these instruments were installed at approximately 2 m above the saplings. Volumetric soil moisture (% by volume) was monitored using Time-Domain Reflectometry (TDR) probes (CS 616, Campbell Scientific, Logan UT, USA), with the default factory calibration, installed in 2 pots for L. speciosa and S. macrophylla and 3 pots for P. indicus.
Tree biometric data. Tree data, including girth at breast height (GBH in cm), and leaf area were also measured to parameterize the model. The GBH was converted to diameter at breast height (DBH), assuming the stems were circular. Sapwood area was estimated from the measured GBH using the equation for calculating a circle, assuming negligible non-conducting wood (pith) due to small stem size (GBH ranged 7.3-10.7 cm). The assumption of negligible pith was confirmed from cutting some stems in the previous experiment and found that 98% of the basal area was sapwood area. For leaf area, we sampled 5 leaves from the upper and lower half of the crown of each sapling, scanned using a scanner (EPSON L3110 model) and analyzed the areas using ImageJ 65 . Total leaf area of each sapling was then determined by multiplying the average measured leaf areas with the number of leaves in the respective crown layer. For scaling purposes using the model, leaf area and sapwood area were determined at the beginning of the study. We realized that leaf shedding may occur in some species, which could change the leaf areas of the studied saplings. However, since the study period was of 6 months, we could not get significant change in leaf area. Also, we considered to preserve the leaves because collecting them for the measurement may influence the sap flux measurements through changing transpirational area. Mean values (± SD) of total leaf area of the saplings were 3 ± 0.34 m 2 for P. indicus, 2.55 ± 0.18 m 2 for L. speciosa and 2.97 ± 0.28 m 2 for S. macrophylla.

Sap flux measurement.
To measure sap flux density (J t ; g m −2 s −1 ), of the saplings, Granier-type thermal dissipation probes 20,21 (TDPs) were constructed. Each T-type thermocouple (copper-constantan) was encapsulated within a steel needle making up the probe. A data logger (used to measure temperature difference between the two probes) was connected to each copper end with the constantan ends connected to each. The sap flux rate www.nature.com/scientificreports/ was measured continuously by tracking the difference between the unheated (upstream or lower probe) and the heated (downstream or upper probe). The temperature difference was converted to sap flux density using the Baseliner program version 4.0 66 . All environmental sensors and sap flux probes were connected to a datalogger (CR1000, Campbell Scientific, Logan, UT, USA) which recorded 30-min average data throughout the study period. However, due to sensor failure for a relatively long time, we could not obtain continuous PAR record, with gaps covering 39% of the period. Consequently, the PAR data were not considered in further analysis after finding only 2% difference between running the model with and without the PAR sub-model.

Model description.
We used the State-Space Canopy Conductance (StaCC) model to gap-fill the missing sap flow measurements using Bayesian inference, and to estimate the latent variables such as canopy stomatal conductance (G t ), transpiration per unit leaf area (E L ), and canopy transpiration (E c ) 31 . The model utilizes a combination of sap flux data and canopy conductance models while accounting for random errors associated with individual sensors. The final hierarchical Bayesian construct is a joint distribution accounting for sap flux measurement, latent states (transpiration and canopy conductance), and model parameters. The data model uses the measured sap flux density to estimate sap flux rate in probe i at time t as: where J t is the average sap flux at time t, Z(d i ) is the sapwood depth sub-model, a i is the random effect due to probe i and S is the observational variance. As the sap flux was measured on potted saplings with small stem size (GBH ~ 7.3-10.7 cm), the sapwood depth sub-model was not used in this study. The model estimates canopy conductance (G t ) and transpiration per unit leaf area (E L ). The steady-state conductance at instant t (G s,t , mmol m −2 s −1 ) is modeled semi-mechanistically 67,68 through multiplicative nonlinear functions of environmental covariates, including vapor pressure deficit (D t ), photosynthetically active radiation (Q t ), and volumetric soil moisture (M t ) as: The sub-model estimating the effect of D t on G s,t has the form: where G ref is the reference conductance (or the canopy conductance at D t = 1 kPa) and λ is the stomatal sensitivity to D t 68 . The dependence of G s,t on light is incorporated through the sub-model: where, α 1 accounts for nighttime conductance and α 2 is the sensitivity to Q t . However, because of the large gaps in PAR record, we did not use this light sub-model, as explained previously. The soil moisture sub-model has the form: where α 3 is the threshold below which M t reduces G s,t and α 4 controls the sensitivity of the reduction in G s,t with decrease in M t below threshold. G s,t is used to calculate actual canopy conductance (G t ) using the state equation: with the assumption that G t depends on the conductance at previous time step t−dt and dt = 30 min and V t = 1 − exp(−dt/τ). The V t term accounts for stomatal lags and τ = 10 min 31 . Canopy conductance is then scaled to calculate E L (kg m −2 s −1 ) as: where q t is a composite variable. We kept a weak prior on the species-specific lags in the capacitance sub-model described in the previous study 31 . Additionally, wherever possible, the same priors for the data and process models were used 31 . A separate model analysis was implemented for each of the three species and the two seasons. The Gibbs sampler was run for 10,000 iterations and the first 3000 iterations were discarded as burn-in. The model generated J t , G t , and E L estimates at each time step and was used to evaluate the differences in water flux among species and across the wet and dry seasons. Daily, weekly, and monthly averages of J t , G t , and E L , were calculated for each species to evaluate the variability in water flux across temporal scales. We then estimated the relationship between daily G t and E L and environmental variables (VPD and soil moisture) using linear regression analyses. The analyses and visualization were carried out on MATLAB 7.12.0 R2017a (The MathWorks, Inc., Natick, Massachusetts, USA). To determine the interspecific and inter-seasonal differences, ANOVA and multiple regressions associated with mean parameters were performed, using the MATLAB program.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.  www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.