Root system variation of pulse crops at early vegetative stage

Grain legumes known as “Pulse crops” are economically and nutritionally one of the most important crop families after cereals. Climate change and abiotic stress factors are limiting yield potential in these crops. Root system architecture, a neglected field, is promising for improved environmental adaptability and higher growth potential. Genotypes with deep and dense root system may cope better with water and nutrient limitations. This study aimed to evaluate 12 common cultivars from six different legume species. Root system architectures were evaluated under controlled conditions in a plexiglass system at the early vegetative stage. Roots were divided into four depth zones (0, 5, 10, 15> cm), and inter and intra species diversity were analysed. Significant diversity was obtained within and between the species. Bean, chickpea, and broad bean constituted deeper and dense root systems while lentil, soybean, and pea formed non-dense and shallower root systems. There was a significant correlation between earliness and early vegetative root vigour. The results of the study may provide a better understanding of the root system architectures of each species-genera. The results presented here may shed light on the selection of root traits in legume breeding programs.


Introduction Introduction Introduction Introduction
The legume family (Fabaceae) contains about 730 genus and 19.400 species and is grown almost worldwide (Rahman and Parvin, 2014). Grain legumes are one of the most important genera for human and animal nutrition and most cultivated product in all but cereals on the world. The most commonly produced grain legumes are dry bean, followed by pea, chickpea, and lentils (FAOSTAT, 2019). They have a vital role due to rich protein composition, and mostly consumed in 3 rd world countries (underdeveloped as well as developing countries), especially in sub-Saharan Africa (Marmot and Bell, 2018), in which nearly 220 million people suffer health problems due to malnutrition (Nedumaran et al., 2015). Legumes are an important part of agriculture as a rotation plant. Their symbiotic interactions with Rhizobium species and nitrogen fixation ability makes them good soil recovery crops (Yahaya et al., 2019). It is estimated that about 140-200 million tonnes of nitrogen are fixed to the soil by biological N fixation and 25 to 30% (35 to 44 million tonnes) of it is fixed by symbiosis annually (Rascio and La Rocca, 2008). Also, it is known that legumes deeply process the soil by their taproot and enable increased aeration (Burridge et al., 2016.). Using legumes in rotation with cereals improves the nutritional potential of soil, prevents soil erosion, and leads to a decrease in the impact of soil pathogens (Daryanto et al., 2015). Microbial activity substantially increases in the rhizosphere of legumes compared to other crops (Lynch, 1990;Houlden et al., 2008).
The two main functions of the root are anchoring, and soil exploration for water and mineral nutrients (Lynch, 1995). Grain legumes need a well-developed root system in the early vegetative stage for high adaptability to stress factors and seed yield (Nicola, 1998). An effective root system can support plants under water and/or nutrient-limited conditions (Chen et al., 2017;Sulis et al., 2019). Some root characteristics of bean and pea species were investigated as early as Weaver and Bruner (1927). Recently, many studies suggested modification of root system architecture (such as deep-root, large root surface area, more effective lateral root distribution), to improve yield, drought tolerance, and tolerance to lack of nutrients (Beebe et al., 2006;Ghanem et al., 2011).
The root architecture (RA) was explained by Lynch (1995) to refer to distinct aspects of the shape of root systems. Ye et al. (2018) also described RA as a significant trait referred by the shape and physical space of the root. RA examines the root system in terms of morphology, topology, distribution, and architecture. The spatial configuration of the root system can be understood from the RA. While root topology examines the root axes, distribution refers to the presence of roots in a positional grid or gradient, each alone is not enough to know the root system. However, root topology and distribution are not needed to know if the RA is known.
Because root architecture offers a 3-dimensional viewing angle (Morris et al., 2017). Generally, soil resources are unevenly distributed or localized in deep soil zones. The RA is the main factor that determines the uptake efficiency of water and nutrients from the soil (Colombi et al., 2018). Moreover, RA refers to competition with the neighbour plant roots and the location of lateral roots.
Climate change and global warming have been speeding up all over the World and various stress factors, especially drought and high temperatures, threatening plant productivity (Santra et al., 2014;Trnka et al., 2015). Novel approaches are needed to cope with changing environmental conditions and to enhance plants' ability to adapt. Various root traits are targeted for adaptation to climate change (Forde, 2009;Ehdaie et al., 2012) such as root tip diameter (Haling et al., 2013;Colombi et al., 2017), gravitropism (Alahmad et al., 2019), and rhizosheaths (Delhaize et al., 2015;Basirat et al., 2019). There is limited knowledge of plants' response to limited conditions and the role of the root system architecture. It is thought that RA can enable us to understand the relationships between some agronomic traits and root system structure. Deeper rooting is a function of both genotype (G) and environment (E), however, the interaction among them makes it difficult to separate genotypic differences from the environmental factors in the RA, due to high G*E interactions and plasticity in roots (Ehdaie et al., 2012). There is a need for understanding the rooting depth and RA under controlled conditions. The outputs of controlled studies would help to identify the right RA forms for the changing climate and soil. Therefore, the aim of this study was to I) evaluate the root system architecture of six different grain legume species, II) to determine inter-and intra-species genotypic variations in early vegetative root traits and III) to investigate the natural growth characteristics of each species under controlled conditions with a modified growth pouch technique.

Used materials and experimental design
The study aimed to evaluate the diversity of early vegetative root system architecture in legume species. Twelve accessions from 6 species (2 cultivars per genus) were used in the study. The used two cultivars in each species were: lentil (Lens culinaris L. cvs. 'Fırat-87' and 'Kafkas'), chickpea (Cicer arietinum L. cvs. 'Arda' and 'Diyar-95'), soybean (Glycine max L. cvs. 'Altınay' and 'Cinsoy'), bean (Phaseolus vulgaris L. cvs. 'Cihan' and 'Yunus-90'), broad bean (Vicia faba L. cvs. , and pea (Pisum sativum L. cvs. 'Green pearl' and 'Sprinter'). The names and some characteristic traits of used cultivars are given in Table 1.  Okçu et al., 2005 The seeds were provided by GAP International Agricultural Research and Training Center (GAPUTAEM), Aegean Agricultural Research Institute, and by some seed trade corporations. The seeds were surface sterilized with 70% ethyl alcohol and 5% sodium hypochlorite (NaClO). After sterilization, the seeds were imbibed in distilled water for 24 hours to initiate homogeneous germination. The next day, seeds were placed on a germination paper (40*40 cm) covered with two layers of plexiglass plates and laid in a container with appr. 45° angle as described in Hohn and Bektas (2020). The plastic container was filled with distilled water up to 70% capacity to provide adequate moisture to seeds. Containers were kept in a climate chamber in fall 2019 for 20 days. The experiment was designed according to randomized complete blocks (replications were treated as blocks) design with four replications and two plants in each replication. The temperature and humidity levels were between 25-27 °C and 60-70%, respectively. On the 20 th day, plates were laid down, and root architecture parameters were measured by getting images of each plant.

Root scanning and image analysis
Root samples were scanned in colored scale at 600 dpi resolution using a portable scanner (ISCAN, handheld scanner). Sample images belonging to scanned plants were given in Figure 1. Each plant is scanned separately and mean values of 8 plants (4 replications x 2 plants per replication) were used for statistical analysis.

Measurements
Root images were analyzed using ImageJ image analysis software (Schneider et al., 2012) to evaluate the phenotypic variability of root architecture. A total of twelve traits constituting root architecture were measured/calculated by the mean values of replications. Root zones were divided into four sub-sections (from the bottom of the stem to root tip). Sections 1, 2, 3 and 4 were representing 0-5, 5-10, 10-15 and 15+ cm depths, respectively. All traits were measured manually for high precision and to capture the value of each root and section (Table 2).    The values calculated were subjected to statistical analysis for the evaluation of significance within and between species. A total of twelve traits (NOR, TLRL, TapRD, aTLRL, NR1, RL1, NR2, RL2, NR3, RL3, NR4, and RL4) were analyzed ( Table 2). The normality of the data was tested using the Shapiro Wilks (1965) normality test. Data transformation was conducted for the data do not fit normality assumptions. Data transformed with logarithmic transformation and fit the normality were evaluated with analysis of variance (ANOVA), while some data did not fit normality were analyzed with the Kruskal Wallis test (Table 4 and 5). Genotypes from the same species were evaluated for the significance within the species. Data with normal distribution within each species were analyzed using a t-test, while the rest of the not fitting normal distribution were analyzed with the Mann Whitney U test (Table 3). One Way ANOVA and Duncan Multivariate tests were calculated using SPSS 21 statistical software (IBM SPSS Statistics for Windows, Version 21.0. Armonk, NY). R statistical analysis software (R ver. 3.5.2, R Core Team, 2018) was used for the Kruskal Wallis test, and the Least Significant Difference (LSD) test was applied to calculate multiple comparison values (De Mendiburu, 2015).

Results Results Results
Twelve genotypes from six different grain legume species were evaluated in a plexiglass early vegetative root analysis system (Hohn and Bektas, 2020). The system allowed us to observe complete roots system architecture and diversity of lateral root initiation within and between legume species (Tables 3-7). Also, the phenotypic observation with this method allowed the measurement of substantial root architectural traits to be collected in an effective, low-cost, and high-throughput design. The plate system, a modified version of the growth pouches, allows identifying main RA traits, with the limitation of mature plant traits under field conditions. There is a need for further studies, looking at the correlation between the early vegetative results and field observations. The characteristics of root system architecture at the early vegetative stage showed considerable diversity at inter and intraspecies levels. The study may help to reveal some unidentified issues about root characteristics and variation in grain legume crops.
According to the results of t and Mann Whitney U tests, statistically significant differences (p=0.05 or 0.01) were determined among cultivars belonging to the same species except for the genus Cicer (Table 3).
There were not any significant differences between the traits of 'Arda' (CA) and 'Diyar-95' (CD) in chickpea cultivars. Significant differences (0.01) were observed concerning all tested traits among 'Kıtık' (VK) and 'Filiz-99' (VF) in broad bean cultivars. The VK and VF exhibited completely different root characteristics. The TapRDs of 'Green pearl' (PIG) and 'Sprinter' (PIS) of the pea genus were significantly (0.05) different than each other (Table 1), however, the rest of the traits were not statistically significant between these cultivars. While there were significant differences in terms of some characteristics, including NOR, LRL, RL2, RL2, NR3, RL3, and NR4, among the varieties in the bean genus, there was no statistical difference between rest of the traits. Although 'Cihan' (PC) and 'Yunus-90' (PY) are in the same genus, they have noteworthy differences in some traits. In terms of aTLRL, RL1, and RL2 traits, significant differences were determined between 'Altınay' (GA) and 'Cinsoy' (GC) cultivars, belonging to the soybean genus, but the other traits were not statistically different in this genus. The lentil cultivars, 'Fırat-87' (LF) and 'Kafkas' (LK) were statistically different in terms of aTLRL, and NR1 traits while others were not. While making a comparison between species-genera, one-way ANOVA was used in the analysis of the traits that showed normal distribution according to the normality test (Table 4), and the Kruskal Wallis test was used for the ones that did not show a normal distribution (Table 5). According to the results of the oneway ANOVA, statistically significant differences (0.01) were determined at NOR, NR1, RL1, TapRD, TLRL, and aTLRL traits. Table 3. Table 3. Table 3.   The significant differences between traits that did not exhibit normal distribution were analyzed by the Kruskal Wallis test (Table 5). In terms of NR2, RL2, NR3, and RL3, statistically significant (0.05) differences were detected while significant differences at p<0.01 level were observed at NR4 and RL4. The highest NOR was obtained on PY (141.7 plant -1 ) while the lowest was on LF (3.7). The PY had the longest TLRL (222 cm) and GC had the shortest one (11.4 cm) among cultivars. The highest aTLRL was obtained on CD (12.6 cm) while the lowest was in GC (0.8 cm). When the number of root distribution in sections was evaluated, the highest NR1 in the S1 section (0-5 cm) was obtained in PY (32.7 plant -1 ) while the lowest value was in LK (5.0 plant -1 ). The longest RL1 measured in the CD (192.8 cm) whereas the shortest one (11.4 cm) was in the GC. The mean values for NOR, NR1, RL1, TapRD, TLRL, and aTLRL are given in Table  6. The traits given below (Table 7) were analyzed by non-parametric methods because the data showed abnormal distribution. The highest NR2 (27.7 plant -1 ) was obtained in the PY while the lowest one (1.0) was determined in the LF. Cultivars in PY and VF showed high potential for the number of roots in the 2 nd section (NR2) and there was no significant difference between them. All cultivars formed lateral roots in section 2, Except the GC. The longest RL2 (116.3 cm) was measured in the VF, followed by PY (31.9 cm). The shortest RL2 (1.5 cm) was recorded in the CD. The 7 cultivars, including GA, LK, PC, PY, PIS, VF, and VK, could constitute lateral roots whereas the rest did not have any root in the 3 rd section. The highest NR3 (27.7 plant -1 ) was observed in the PY while the lowest one (0.3) was in the LK. As the length of secondary roots in the 3 rd section, the longest RL3 (34.9 cm) was measured in the VF while the shortest (0.4 cm) one was determined in the LK. Although half of the cultivars had secondary roots in the 4 th section, only two cultivars had a secondary root that is longer than 2 cm. So, it was recorded that most cultivars do not have a noteworthy secondary root after the 3 rd section in the early vegetative stage. The highest NR4 (53.7 plant -1 ) was determined in the PY whereas the lowest value (1.3) was obtained from the PC. The longest and shortest RL4 (33.0 and 0.2 cm) were measured in the PY and GA, respectively (Table 7).
All species constituted highly numerous lateral roots except lentil cultivars in section 1. The highest number of lateral roots were observed in the bean genus, followed by the VF and GA. The root numbers exhibited significant differences between genera and even species. This event is an indicator that the early vegetative stage root growth rates of the species are different from each other and that the root growth rate in lentil is slower than other species. In contrast, common bean species exhibited the highest growth rate in not only the 1 st section and also the other parts of the root.  In 20 days, lateral root formation of species was intensely observed in section 1 and it strongly decreased in other sections in all but the PY and VF. Lateral root formation continued to increase up to section 4 in the PY and achieved the top level. It was determined that cultivars either started to form lateral roots yet or failed to form any roots in the 4 th section, except the PY and VF. Another point to note is that the growth rate did not only differ between species but also exhibited considerable differences among cultivars belonging to the same species (Figure 2). It was seen that seed size (endosperm size) had a significant effect on early vegetative root development, as expected. This effect is confirmed with six different legume species in our study.  The root length density was more effectively observed in section 1, as well as the number of roots. However, although number of roots was determined low in some cultivars, the total root length for section 1 was observed higher. The longest RL1 was obtained on the PY, followed by the VF and PC. Although there were nearly no differences in RL1 values, except PIG and PIS, statistical differences were determined in terms of NR1 among cultivars that are from the same species. From a general perspective, the highest root length values were determined in broad bean and bean while the lowest ones were obtained in lentil and soybean. The longest RL2 was determined in the PY and VF whereas the shortest one was seen in the LF. Also, the GC did not form any root in section 2 (Figure2). A longer total root length in VF was measured compared to the PY in section 2, although NR2 of PY was higher compared to VF. This is an indicator that the PY has a higher root formation capacity, but the root growth rate of the VF was superior in section 2. The RL value of the CD, which had the longest RL1, strongly decreased in section 2. The root length of the CD was gradually decreased in deeper layers (Figure 3).

Discussion Discussion Discussion
The experiment was conducted in controlled conditions using a modified plexiglass system. The root characteristics of 12 cultivars were investigated and a wide genetic variation was observed within and between species (Table 3, 4, and 5). Non-destructive and easy screening of the RA, easy root development observation, time-saving processing, easy repeatability, and no need for any root washing is possible with the plexiglass system. However, the used technique can not completely reflect field conditions due to eliminating interaction with beneficial microorganisms. The controlled conditions also inhibit plant's exposure to environmental factors; therefore, the physiological relevance of roots needs further evaluation under field conditions (Pierret et al., 2013;Wasaya et al., 2018).
The study showed that nearly every genotype exhibited different root structures. A wide variation was determined within and between six species. The RA is determined and controlled by some substantial traits, including the number of total and lateral roots, total root length, taproot length, root biomass, and sections of roots. Also, it was stated that the root system architecture is especially determined by the number of roots, length of the taproot, and length of lateral roots (Mia et al., 1996).
Cultivar level comparisons within each species, yielded significant differences on the broad bean, while cultivars in some species such as chickpea did not show any significant difference. This might be due to the low genetic variation and pedigree similarities within the same species. These cultivars, generally selected for specific target regions, if both cultivars were selected/breed for similar environmental conditions, they would have similar agronomic properties. Similarly, we did not observe any significant differences within lentil and pea. These cultivars in each species exhibited substantially similar characteristics. However, the large differences among broad bean varieties are thought to be affected by the contrasting agronomic properties of the two cultivars (Table 3). Many studies showed that intra-species genetic variation increases with an increasing number of used genotypes. Some studies reported a wide intra-species variation in lentil (Idrissi et al., 2015), bean (Teran et al., 2019), pea (Mcphee, 2005), broad bean (Belachew et al., 2018), chickpea (Chen et al., 2017 and soybean (Suematsu et al., 2017). Therefore, the reason for the low intra-species variation observed in some species may be due to the low number of varieties used.
The longest TapRD values were observed in bean and broad bean while the shortest ones were determined in lentil and soybean (Table 6). This difference might give us an opinion about the early-stage root growth and resource requirement of species. Different root structures provide some advantages in water and nutrient uptake, tolerance to various environmental stress factors, and adaptation to different regions. Evaluation of the root system diversity and architecture in crop genotypes is a pre-requisite for the use in breeding (Paez-Garcia et al., 2015;Kuijken et al., 2015). Taproots give anchorage to plants in legumes and constitute the RSA, control rooting depth, and determine the ability to grow under compacted soil layers. Taproot length is closely related to tolerance in drought stress. A deep root system provides an advantage to benefit from the water in the deep soil layers when there is water deficiency. While lateral roots are the main organ in water and nutrient uptake under optimum conditions, taproot plays a more effective role in drought conditions (Lindh et al., 2014). Sponchiado et al. (1989) revealed that the plants with a deeper root system have higher grain yield and are more tolerant of drought under water deficiency. White and Castillo (1988) found that drought-sensitive bean genotypes could not reach even 0.8 m under drought while the rooting depth of drought-tolerant bean genotypes could reach up to 1.2 m. A deep root system can supply some advantages to soybean (Mwenye et al., 2016), chickpea (Varshney et al., 2011, lentil (Rabani, 2018), and pea (Couchoud et al., 2020) under drought conditions.
The NOR and TLRL allow shaping the rooting potential of plants (Canadell et al., 1996;Schenk and Jackson, 2002). The relationship between NOR and TLRL is still very not clearly defined. Although total root length increases with a higher number of roots it is not enough alone. The maximum NOR value (141.7) was determined in PY while the longest TLRL (361.8 cm) was measured in VF (Table 6). So, to achieve maximum root length both high NOR and long TLRL may be needed. Lateral roots are more effective underground organs due to its role in surface water and mineral uptake that may constitute a substantial part of the RA. The primary targets of lateral roots, especially fine roots, are to extract water and nutrients for plants (Rostamza et al., 2013). Mia et al. (1996) reported the highest total root length in chickpea, followed by cowpea and grass pea in their study with six grain legume species. They also stated that mung bean, lentil, and black gram exhibited a shorter root system. While the bean and broad bean varieties exhibited superior characteristics in terms of NOR and TapRD, pea and lentil showed a shorter root system in our study. The ability of roots to effectively and quickly occupy the soil volume during the early growth stages can be a vital trait related to soil water and nutrient uptake efficiency (Jamont et al., 2013). A well-distributed root system in the soil allows for higher N uptake and it provides some advantages to plant when competing with other plants/species (Dayoub et al., 2017). Thus, a deep root system with lateral expansion might be a key factor affecting the competitive ability of legumes for soil resources during the early growth stage. Furthermore, species that have a greater early growth ability are generally more competitive for soil resources in their life cycle (Andersen et al., 2007). The reason for rapid lateral root formation, especially in erect growing and species with high stature, can be associated with the desire of the plant to anchor onto the soil. Also, compared to slow-growing species such as lentils, the need for water and nutrient uptake supports this hypothesis. The NR exhibited a noteworthy distribution among cultivars in the present study. Root density generally was more shown in the 1 st section followed by the 2 nd , 3 rd, and 4 th and showed similarity with RL. Especially the PY and VF came into prominence in terms of root distribution at deeper sections. By the 20 th day, some accessions did not have yet formed lateral roots in deeper sections (Section 3-4), or newly formed (Figures 2 and 3). The root length density on different sections (depths) of the taproot is a substantial trait that is one of the most important factors affecting water and mineral uptake. Because different mineral nutrients disperse in particular soil layers and the restriction rate of them showed differences for plants. It was stated that the concentration of some nutrients, especially N, P, and K, are higher than the others such as Na and Cl that usually restricted to the upper 20 cm of soil (Jabbagy and Jackson, 2001). So, increasing lateral roots provides an advantage to plant to benefit from the macro-nutrients. The number and length of roots on different sections of taproot were shown in Figures 2 and 3, respectively. Distribution and density of roots in soil layers change depending on genetic composition, growth stage, soil structure, and environmental factors/limitations (Manschadi et al., 2013). The water uptake pattern is uneven throughout the depth of the rooting zone. The water is taken by plants changes depending on soil structure, the existence of water in the soil, root architecture, root density, and different zones of taproot (Pardo et al., 2000;Wasson et al., 2012;Zhang et al., 2020). In general, almost 40% of the total water uptake occurs over the 1 st one-fourth of the root zone, 30% over the 2 nd , 20% over the 3 rd , and 10% over the last 4 th of the total rooting depth under optimum water conditions (Subbarao et al., 1995). Water and mineral uptake ability from different soil depths provides an advantage especially when there is a resource limitation. However, plants often tend to form higher root length densities in the surface layer than in the deeper zones. The main reason for this situation can be summarized as follow: I) Quick uptake of water before evaporation after precipitation, II) Creating a reserve against problems caused by disease or pests, III) Reaching inactive nutrients such as P, which are more prominent in the upper layers of soil, and IV) To provide and advantage in the competition of minerals on surface layers (Passioura, 1983;Subbarao et al., 1995).
A significant difference in the number of roots in each section was seen in species and cultivar level. Detection of a wide variation among varieties belonging to the same species expresses the diversity of root traits. Weaver and Bruner (1927) stated shallower root development in pea, and little or no lateral rooting at depth. Our results are correlated with Weavers' observations. Two broad bean cultivars had quite different rooting patterns, VF had deep rooting with many laterals, while VK had poor rooting ability and almost no laterals. When the agronomic properties of the broad bean cultivars are compared, VF has earliness, while VK matures later (Table 1). It was seen that a cultivar with the earliness seems to have much faster root initiation. The same earliness-rooting ability interaction was also seen among the bean cultivars. The PY with earliness characteristic seems to develop a more advanced root system in an earlier stage than PC, which did not have any meaningful lateral roots in the 4 th section (15+ cm depth). Additionally, both PY and VF had more lateral roots in the 4 th section than in the 3 rd section. Pereira et al. (2016) stated that peanut cultivars with earliness had 19% more root length under water-stressed conditions and also had a higher rate of root formation under water stress compared to cultivars with longer maturity duration. On the other hand, Subbarao et al. (1995) suggested that pigeonpea cultivars with later maturing characteristics developed deeper roots compared to early maturing ones. However, these results were reported from physiological maturity. It is a predictable hypothesis that early maturing varieties tend to form a root system more quickly to adapt to the environment and gain an advantage in competition, but those late-maturing varieties may have the potential to create a more effective root system eventually. To better understand the relationships between early root formation and genetic factors, there is a need for more detailed studies in this field. It is thought that a root system, which is rapidly established in deeper soil zones, will provide an advantage in competition with neighboring plants in terms of utilizing resources in different soil layers. Distribution of lateral roots in different depths of the taproot is a noteworthy character to determine the rhizosphere coverage area and to exploit the water and nutrient composition of the soil (Lauenroth and Gill, 2003).

Conclusions Conclusions Conclusions Conclusions
The study showed that there is a wide variation both at the genus and cultivar levels in grain legume species. We were able to observe the genetic potential of the species and cultivars without the limitation of the genotype*environment interactions. Each species demonstrated unique root architectures under controlled conditions. We observed a correlation between earliness and rooting ability. Cultivars with earliness tend to have more advanced roots (higher number of laterals and longer total root length) at the early vegetative stage. Understanding earliness and root development interactions would allow breeders to make an indirect selection for roots when breeding cultivars for drought and/or nutrient-limited conditions. However, there is a need for further evaluation of various grain legume cultivars under field conditions to interact with these results with field observations.

Authors' Contributions Authors' Contributions Authors' Contributions Authors' Contributions
MC and HB designed, run the experiments, collected data, and wrote Ms draft; FY made statistical analysis and contributed writing, ME edited Ms and gave feedback. All authors read and approved the final manuscript.