Harvest effects on density and biomass of Neopicrorhiza scrophulariiflora vary along environmental gradients in the Nepalese Himalayas

Abstract A surprisingly large number of species potentially threatened by human harvest lack quantitative ecological studies incorporating harvest effects, especially clonal species in the alpine Himalayas. We studied density and biomass variation of a threatened medicinal herb, Neopicrorhiza scrophulariiflora, to examine the effect of harvest on plant performance. The study covered two regions with contrasting harvest situations—one with open‐access and another protected from commercial harvesting. Four populations from each region were compared along an elevation gradient (3,800–4,800 m). Also, we conducted in situ interviews with 165 and 38 medicinal and aromatic plant users in open‐access and protected regions, respectively, to assess the collection and use patterns of the target species. The quantity harvested per household for traditional healthcare use was similar in both regions. We found no evidence of trade‐driven collection in the protected region but in the open‐access region a trade‐based annual collection of 35–465 kg dried rhizomes per household had a strong negative effect on both density and biomass. In the protected region, the effect of harvest intensity on plant density was positive for vegetative and negative for reproductive individuals, whereas in the open‐access region, the effect was negative for both vegetative and reproductive individuals. The results indicated that a low harvest intensity had no adverse impact on N. scrophulariiflora populations; however, quantification of the optimum level of harvest remains to be explored. Shrub vegetation appeared to buffer the harvest impact on plant density, possibly through the retention of additional moisture. To maintain population viability, we suggest regulating harvest, for example, by introducing rotational harvest systems, ensuring that a sufficient number of reproductive individuals are left as a source of propagules in each harvested population and that populations are given time to recover between harvests.

and physical damage are important services that can be offered by the surrounding vegetation (Ballantyne & Pickering, 2015;Callaway et al., 2002). Thus, the balance between facilitative and competitive environmental interactions has significant ecological implications to the performance of alpine plant species (Chu et al., 2009;Klimešová et al., 2012).
In addition to environmental interactions, alpine plants are in many cases also subjected to high levels of anthropogenic impacts. In the Himalayas, for example, international trade in medicinal and aromatic plants (MAPs) (Olsen, 2005;Pyakurel, Sharma, & Smith-Hall, 2018) has created extreme pressure on populations of some alpine plants and their habitats. In addition, other factors, such as livestock grazing and fire, may reduce densities of alpine plant populations (Niu, He, Zhang, & Lechowicz, 2016). Harvesting of whole plants or plant parts may affect reproduction, survival, and growth, thereby affecting plant population dynamics (Gaoue, Horvitz, Ticktin, Steiner, & Tuljapurkar, 2013;Ghimire et al., 2008;Ghimire, McKey, & Aumeeruddy-Thomas, 2005;Huai, Wen, Xu, & Bai, 2013;Ticktin, 2004). The extent of harvest impact on plant populations, however, varies depending on habitat conditions, plant growth strategies, and species-specific processes, such as regeneration and colonization (Gaoue et al., 2013;Ticktin, 2004Ticktin, , 2015. Clonal plant species to some extent buffer the anthropogenic effects like harvesting by mobilizing stored reserves to maintain normal phases of growth (Mandle & Ticktin, 2012). Furthermore, the clonal tendency toward developing a bud bank after a disturbance also plays a crucial role in vegetative reproduction (Evette, Bédécarrats, & Bornette, 2009).
The clonal architecture, especially in plants exhibiting a "guerilla" strategy of clonal growth, producing long and spreading vegetative offshoots, allows the plant to escape stressful microsites and colonize more suitable ones (Ghimire et al., 2005;Humphrey & Pyke, 1998). In clonal herbs exhibiting a "guerilla" strategy, the density of plants is positively affected by light harvesting if it occurs at sufficiently long intervals and does not exceed the regeneration potential (Ghimire et al., 2005). Thus, with responsible management practices, vegetative propagation in clonal plant species could compensate for reduced biomass caused by anthropogenic disturbances.

The interactions between people and plants in the Himalayas
are not yet fully understood (but see Ghimire et al., 2005;Ghimire et al., 2008;Rokaya, Münzbergová, & Dostálek, 2017). Specifically, determining the relative effects of habitat factors, such as substrate, vegetation cover, and disturbance on plant density and biomass is crucial for development of well-informed and effective conservation policies. This study aims to disentangle the factors that control the density and biomass of a clonal medicinal herb, Neopicrorhiza scrophulariiflora (Pennell) D.Y. Hong, across a continuum of environmental conditions and harvest intensities in the alpine Himalayas, Nepal. This species is highly threatened due to overharvesting of its rhizome for regional and international trade (Ghimire et al., 2005). Studies emphasizing the variation in density and biomass of such species across a continuum of environmental conditions along an elevation gradient could enhance our understanding of the ability of the species to respond to environmental stress. We will address three questions: (a) In what ways do plant utilization patterns vary between study regions in terms of the purpose and intensity of harvest, and whether harvest takes place in the beginning or at the end of the growing season? (b) How do elevation, surface and vegetation cover, and anthropogenic disturbance influence the density and biomass of N. scrophulariiflora?
and (c) Does plant response to harvest impact vary between study regions? To answer these questions, we first assessed plant utilization systems and analyzed density, population structure, and biomass variation in different populations along elevation gradients in two study areas characterized by different protection regimes. The species' interaction with environmental factors, including habitat conditions, surrounding vegetation, and human harvest regimes, was analyzed to evaluate these associations.

| Study area
The study was conducted in two protected areas (hereafter re- The climate in both regions varies from subtropical to alpine-nival. The two regions are similar in terms of annual precipitation (average of 2,100 mm in ANCA, and 2,300 mm in LNP) and annual average temperature (range: 4-27°C for ANCA and 2-25°C for LNP; both precipitation and temperature were recorded at the nearest meteorological stations located at 1,900 and 1,100 m above sea level, respectively). Both regions are characterized by high mountains and rugged terrain. Almost 70% of residents in ANCA and about 50% in LNP live below the poverty line Uprety, Poudel, Asselin, Boon, & Shrestha, 2011). Local livelihoods depend on traditional agriculture, livestock rearing, and seasonal labor (CBS, 2012). The upper slopes are fragile and arable fields are limited to low elevations (<3,000 m), but the productivity is lower in ANCA than in LNP (CBS, 2012;DDC, 2015;Fox, Podger, & Yonzon, 1996). Therefore, people in ANCA supplement their income mostly by harvesting and selling high-valued MAPs and other forest products (Pouliot, Pyakurel, & Smith-Hall, 2018;Pyakurel et al., 2018). In LNP, people primarily supplement their income through livestock rearing, cheese production, and tourism, and are less involved in harvesting and trade of forest products (Shakya, 2016 and LNP as "protected." F I G U R E 1 Map of the study area: triangles indicate populations included for the ecological study in Api-Nampa Conservation Area (ANCA, left) and Langtang National Park (LNP, right). The thin and thick contour lines are provided to define the elevation scale in every interval of 500 m and 1,000 m from 3,000 to 5,000 m (source of land cover map: Uddin et al., 2015)

| Study species
Neopicrorhiza scrophulariiflora (Pennell) D.Y. Hong (hereinafter referred to as Neopicrorhiza) of the family Plantaginaceae (previously Scrophulariaceae) (APG IV, 2016) is a perennial, rhizomatous herb that exhibits clonal growth. It is confined to the pan-Himalayas, the Tibetan Plateau, Assam-Burma, and south-central China, and is found at elevations between 3,500 and 4,800 m (Ghimire et al., 2005;Hong, 1998). The plant is relatively slow growing and is mostly confined to cool north-facing slopes (Shrestha & Jha, 2009).
Flowering and fruiting occur in June-September. Seedling recruitment is low, and propagation mainly takes place by vegetative means (Ghimire et al., 2005). The modular growth of the individual genet takes place by horizontal extension of ramets (vegetative offshoots), which can survive independently when detached. The clonal growth of Neopicrorhiza has been described as a "guerrilla" strategy, leading to spaced clusters of spreading fugitive ramets connected below ground by horizontal rhizomes (Ghimire et al., 2005).
The rhizomes of Neopicrorhiza are valued for the treatment of cough, cold, headache, fever, high blood pressure, conjunctivitis, bile reflux, and intestinal pains, among others (Ghimire et al., 2005;Manandhar, 2002). It is one of the most commonly consumed herbs in the Himalayas and is appraised to have a high pharmacological potential (Li, Liu, Abdulla, Aisa, & Suo, 2014). The plant is resistant to grazing, and thus, premature and excessive harvesting is the main issue for its conservation (Ghimire et al., 2005;Shrestha & Jha, 2009).
It is closely related to Picrorhiza kurrooa Royle, a species native to north-west India listed in CITES Appendix II. Rhizomes of both species are traded globally, but Neopicrorhiza is more commonly traded than P. kurrooa (Mulliken, 2000). Rhizomes of Neopicrorhiza collected from Nepal are mostly traded to India in air-dried form through wellestablished marketing chains (Olsen, 2005) with a current annual trade volume of ca. 350 tons (personal communication with traders in ANCA and at a major road-head collection center in Nepalgunj, south-west Nepal, November/December 2017).

| Interviews with local resource users
This work was carried out as a partial basis for the first author's Ph.D. dissertation at the Central Department of Botany (CDB), Tribhuvan University, Nepal. The study proposal was approved by the CDB Research Committee. Written permission for fieldwork was obtained from Department of National Park and Wildlife Conservation, Government of Nepal. The research objectives were explained, and prior informed consent was received from respondents before conducting interviews and consultations with MAP users in ANCA and LNP. The participants were informed they could withdraw their consent at any point in case of disagreement.
Commercial MAP harvesting and livestock grazing are the main human activities in the alpine regions of the Himalayas. Local people usually follow traditional systems of rotational grazing when livestock are brought to subalpine and alpine pastures for summer grazing. We recorded two types of MAP collectors: dedicated collectors, who collect only one target species at a time and are not involved in pastoralism; and opportunistic collectors, who collect different MAP species and are also involved in pastoralism (Olsen, 1998

| Plot establishment and sampling
We carried out a sample-based study in Neopicrorhiza populations dur-  Table   S1). In each population, we laid out three transects at intervals of 100 m (Ghimire et al., 2008(Ghimire et al., , 2005. In each transect, we established two plots (3 × 3 m) at horizontal distances of at least 20 m from each other, depending on field conditions. The clonal behavior of the study species is such that sprouting occurs readily from the parental ramet within a limited distance leading to a patchy form of distribution. Therefore, we established the plots subjectively where the plants occurred in sufficient number for sampling (Whittaker, 1956). We divided each plot into nine subplots (1 × 1 m). Among these, five subplots (four at the corners and one at the center) were systematically sampled and measured. In essence, the study included 48 plots with 240 subplots, and each of the samples in the two regions included half of these.
We studied plant density for ramet stage classes. In each subplot, we classified individual ramets, based on number and size of leaves and occurrence of reproductive organs, into one of the three stage classes (Oostermeijer, Brugman, Boer, & Nijs, 1996;Weppler, Stoll, & Stöcklin, 2006): (a) juvenile ramets bearing 1-4 smaller-sized (range: 0.5-2.0 cm length) leaves; (b) adult vegetative ramets bearing more than four larger-sized (>1.5-9.0 cm length) leaves but with no indication of sexual reproduction; and (c) reproductive ramets bearing any number and size of leaves and a reproductive peduncle. To describe population structure, we calculated the proportion of ramets in each stage class.
For the estimation of harvestable rhizome biomass, we asked local harvesters to collect 2-6 full-grown ramets (adult vegetative and reproductive stages only) from eight to eighteen 1 × 1 m plots in each population, laid out randomly in the vicinity of plots used in density estimation. This way we collected biomass samples from 50 subplots in ANCA and 46 in LNP. In total, we extracted 161 biomass samples in ANCA and 133 in LNP (30-59 samples including above-and below-ground parts per population). We measured the length and girth of each rhizome when fresh. Girth was calculated as a mean of three measurements at different nodal points along the rhizome. The volume of each rhizome was calculated based on girth and length assuming cylindrical shape. The biomass of each rhizome sample was measured after drying for 72 hr at 60°C in a hot air oven.
The volume and biomass were used for calculation of rhizome tissue density (expressed in g/cm 3 ). The total number of rhizomes per kilogram dry weight was calculated for each population based on the estimated mean biomass.

| Plot-based environmental variables
Geographical location (latitude and longitude) and topographic features (elevation, slope, and aspect) were measured for each plot.
The shrub canopy cover was estimated visually as total canopy cover of all shrub species, and shrub height was measured at three points located diagonally across the subplot. Cover of ground vegetation was estimated separately for herbs (all dicotyledons species, orchids, and lilies), graminoids (including all grasses and sedges), and lichens/ mosses. Similarly, bare ground cover and rock cover were estimated for each subplot. Soil pH was measured using a portable pH meter Evidence of MAP harvesting and livestock grazing was confirmed by intensive observation in each subplot. In order to determine the magnitude of impact, subplots were further divided into four cells, each 0.25 m 2 in size. We assigned a binary score to each cell based on whether there was evidence of disturbance (1) or not (0). At the subplot level, the scores given separately for each category of disturbance were combined for the four cells to obtain an overall score ranging from 0 (null) to 4 (very high). The following types of evidence were used in assessing the harvesting of Neopicrorhiza: the presence of left-over rhizome fragments, uprooted plants, and holes or scars left after the excavation of plants (Figure 2). Assessment of grazing was based on the presence of pugmarks and trampling, bite marks, browsed plants, broken, defoliated or disorientated aerial parts and animal droppings.

| Data analysis
Variation in ramet density, biomass, and environmental covariates among populations and between regions was tested by nonparametric Kruskal-Wallis one-way ANOVA and Mann-Whitney U tests. A generalized linear mixed effects model (GLMM) was used to examine the effect of environmental factors on ramet density, and a linear mixed effects model (LMM) was used for analyzing rhizome biomass and tissue density. A negative binomial distribution was assumed for ramet density, and rhizome biomass and tissue density were transformed using the natural logarithm to make their residuals normally distributed. To control for intra-class correlation caused by the nested design, "plot" was included as a random effect in models for ramet density. Similarly, a random effect of "subplot" was included in models for rhizome biomass and tissue density. Effects of region and population were treated as fixed nominal factors. All numeric environmental predictors (Table   S2) were treated as covariates. These covariates were standardized (mean 0, standard deviation 1) to avoid erroneous conclusions when analyzing interactions between continuous variables and to facilitate the interpretation of variables and models (Bråthen & Lortie, 2016). We used qq-plots to evaluate the properties of residuals in models based on different distribution assumptions (Zuur, Ieno, Walker, Saveliev, & Smith, 2009). In the case of ramet density, the primary additive model was prepared by starting with F I G U R E 2 Neopicrorhiza scrophulariiflora growing in (a) heavily harvested and (b) unharvested habitats; an individual plant at fruiting stage is 5-25 cm tall. The approximate size of objects is indicated by the scale bars

| Variation in plant density
Ramet density (m −2 ) varied within and between the regions ( Figure 3). Importantly, for total ramet density (all stages combined),

| Effect of environmental factors on plant density and biomass
The environmental variables recorded in the sample plots varied within and between regions (Table S2). Major differences between regions were observed for ground cover variables and disturbances.
Shrub, rock, and moss/lichen cover were higher in LNP than in ANCA. Cover of herbs and graminoids was higher in ANCA than in LNP. The average scores of harvesting and grazing were remarkably higher in plots studied in ANCA than in LNP (Table S2).
Out of twelve variables (including "population") studied (Table   S2), only nine were used in GLMM after removing confounded variables. Moss/lichen cover was confounded with shrub cover, and bare ground cover was confounded with herb and graminoid cover. Shrub height revealed several missing values and so, it was excluded in the model. Generalized linear mixed effects model identified five variables that exhibited a clear association with plant density (Table 1). Rock cover, in general, did not show any notable effect on plant density when data from both regions were treated together. However, when data were treated separately for each region, rock cover was found to be associated positively with plant density in ANCA (Pearson correlation, r = 0.323, 0.281, 0.320 and 0.311 for juvenile, adult vegetative, reproductive and total ramet density, respectively; p < 0.010; Figure 5a). By contrast, in LNP the corresponding relationship was negative (r = −0.303 and −0.275 for adult vegetative and total ramet density, respectively; p < 0.010; Figure 5b).
In the combined model, shrub cover also appeared to have a positive effect on ramet density (Table 1). When data were treated separately for each region, it emerged that increasing shrub cover was associated with increasing ramet density in LNP (Figure 5d)  Among all environmental variables studied, harvest disturbance had a negative effect on ramet density for all stage classes, whereas the effect of grazing was positive for juvenile and adult vegetative ramet densities (Table 1). A significant interaction between harvest impact and region further indicated that the negative effect of harvest on plant density was stronger in ANCA, and the effect of harvest was subtle for the LNP populations ( Table 2). The weak positive effect of harvest on ramet density in LNP (positive interaction in Table 2) for all vegetative stages could be linked to the narrow range of harvest intensities.
The harvest impact showed different association patterns with rock and shrub cover in the two regions. In LNP, the harvest score  Note: Models were estimated assuming a negative binomial distribution and "plot" was included as a random factor. Levels of significance are indicated by the symbols: p < 0.0001 "***," p < 0.001 "**," and p < 0.05 "*." relationship was positive in ANCA (r = 0.326, p < 0.001). Contrary to this, in ANCA, the harvest impact score was negatively related with rock cover (r = −0.268, p = 0.003); whereas, in LNP, the relation was not significant. The positive relationship between harvest impact score and shrub cover further indicated that shrubland habitats in ANCA were preferred collecting sites for Neopicrorhiza harvesters.
In general, while fitting variables in GLMM (  (Table 2). Elevation generally had negative parameter estimates, indicating that ramet density gradually decreases with increasing elevation. In contrast to the other stages, reproductive ramet density in LNP was weakly positively related to elevation (Table 2) (Table 2). We did not observe an effect of elevation on biomass (Table 2), but rhizome tissue density increased with elevation in both regions (Table 2). This indicates that tissue density was more sensitive to environmental variation along the elevation gradient than biomass. Additionally, this may suggest that older ramets with F I G U R E 5 Relationship between total ramet density (all stages, m −2 ) and (a, b) rock cover, and (c, d) shrub cover for Neopicrorhiza scrophulariiflora in Api-Nampa Conservation Area (  Note: Parameters were estimated assuming a negative binomial distribution for density, and natural log transformation was used to control for the skewed distribution of biomass and tissue density. A random effect of "Plot" was included for ramet density, and for biomass and tissue density of rhizome "Subplot" was included as a random factor in the models. Levels of significance are indicated by the symbols: p < 0.0001 "***," p < 0.001 "**," and p < 0.05 "*." denser tissues were more prevalent at higher elevation, especially in populations in LNP (Figure 6a,b).

| Utilization practices and pressure on wild MAP resources
We observed less harvest pressure in LNP (which is protected) In Nepal, Neopicrorhiza is one of the high-valued MAPs of the alpine region with a century-long trade history (Olsen, 2005;Olsen & Larsen, 2003). In recent years, the exploitation of herbal and fungal MAPs has become a significant challenge for sustainable resource management in the Himalayas (Cunningham et al., 2018;Negi, Pant, Joshi, & Bohra, 2014). Presently, due to lucrative markets and surging prices, the collection of the highly valued parasitic fungus Ophiocordyceps has exceeded the total income from the harvest of all other MAPs (Negi et al., 2014;Pouliot et al., 2018). However, when valuable MAPs become scarce, collectors start collecting other highly valued products (Lange, 2004;Mulliken, 2000).
Since livelihood options in ANCA are constrained by its rugged terrain and remoteness, people supplement their income primarily by harvesting MAPs and other forest products (Pouliot et al., 2018;Pyakurel et al., 2018). In contrast, LNP is an attractive tourist destination in Nepal, in part because it features an important pilgrimage site (the sacred Gosainkunda lake), which is visited by tens of thousands of devotees every year (Shakya, 2016). People in LNP therefore supplement their income by engaging in tourism and hotel business instead of harvesting and trading forest and environmental products.
The broad geographic distribution of Neopicrorhiza in the Himalayas implies that this species is an easily accessible economic resource to many rural people (Olsen, 2005;Olsen & Larsen, 2003). The

cash income earned by selling high-value medicinal plants, including
Neopicrorhiza can be an appreciable supplement to local subsistence (Olsen, 2005 Nevertheless, some studies indicate that illegal harvest of MAPs could occur in both regions due to geographical remoteness and lack of regular invigilation mechanisms Shrestha, Prasai, Shrestha, Shrestha, & Zhang, 2014;Uprety et al., 2011).

| Variation in plant density and biomass
Our observations show that the variations in ramet density and rhizome biomass of Neopicrorhiza correlate with harvest F I G U R E 6 Relationship between elevation and rhizome tissue density of Neopicrorhiza scrophulariiflora in (a) Api-Nampa Conservation Area (ANCA) and (

| Major habitat determinants
In LNP, the shrub cover and moss/lichen cover showed a positive relationship with Neopicrorhiza density (r = 0.787, p < 0.0001).
Hence, it could be suggested that shrub vegetation with substantial moss/lichen ground cover maintains a high moisture level that favors Neopicrorhiza populations and leads to increased density.
Positive effects of shrub facilitation, such as reduced evapotranspiration, improved thermoregulation, and soil amendment by litter accumulation and decomposition have been reported in previous studies (e.g., Ballantyne & Pickering, 2015). However, in LNP the shrub cover generally decreased with increasing elevation. Based on the density and biomass estimates in LNP, we suggest low elevation populations (<4,500 m) are likely more resilient to low-intensity harvest, in part due to facilitative effects offered by shrub cover.
However, in ANCA, rock cover showed a positive relationship with plant density (Figure 5a), whereas the relationship between shrub cover and plant density was not significant. Our interviews with MAP collectors revealed that their impression was the opposite, as they believed that, compared to other habitats, shrubland habitats had higher biomass of Neopicrorhiza per unit area, which allowed them to collect larger amounts within a given period of time (data not shown).
Commercial collectors usually harvest intensively in order to maximize yield while minimizing effort and time. Hence, they prefer to harvest in dense populations where high quantities of rhizomes are available (Ghimire, McKey, & Aumeeruddy-Thomas, 2004). We therefore suggest that the observed higher density of Neopicrorhiza in rocky habitats in ANCA may be due to the preference of harvesters for the shrubland habitats at lower elevations.
In our study, grazing showed a positive effect on ramet density except for the reproductive stage (Table 1). Our field observations indicated that Neopicrorhiza is unpalatable to herbivores, which could be linked to secondary metabolites, such as terpenoids, cucurbitacin, kutkisterol, and steroids present in the plant (Li et al., 2014).
Additionally, in alpine meadows, grazing tends to decrease competition for light without altering functional diversity of foliar traits (Niu et al., 2016). Hence, the removal of species other than Neopicrorhiza by grazing reduces competition for light and provides the species with suitable conditions for vegetative regeneration. Indeed, the guerrilla growth form of Neopicrorhiza allows it to multiply as soon as it finds open space (Ghimire et al., 2005;Humphrey & Pyke, 1998).
Based on the parameter estimates obtained for juvenile and adult vegetative ramet densities it appears that for given levels of harvest and elevation, densities were higher in LNP than ANCA (Table 2). Ghimire et al. (2005) also found higher values of Neopicrorhiza density at low levels of harvesting. The plagiotropic growth of the plant enhances its ability to adapt to extreme environmental conditions where the competition from orthotropic plants is less important (Ghimire et al., 2005;Katayama et al., 2015). The ramets are only partly autonomous, which enables them to benefit from nutrients (nitrogen, phosphorus, potash, and urea) and photosynthates (sugars, amino acids, and hormones) made available by already existing ramets of the same genet (Lei, 2010). Previous studies (Evette et al., 2009;Johansen, Wehn, & Hovstad, 2016;Lei, 2010) found that the budding potential and ramet production of clonal herbs in disturbed environments is higher than in undisturbed environments. Consequently, the clonal growth and fugitive establishment of Neopicrorhiza could partly remediate the impact of low-intensity harvest in disturbed populations.

| CON CLUS ION
The present study provides baseline evidence on plant utilization systems, population differentiation, and habitat conditions of

ACK N OWLED G M ENTS
We would like to thank Santosh Thapa, Bikram Jnawali, and Amar Bista for their assistance during fieldwork. We are obliged to Dipesh Pyakurel and James Lucas for their constructive comments on an earlier version of the manuscript. We are thankful to local communities of LNP and ANCA for their support and for sharing their knowledge.
We acknowledge the Department of National Parks and Wildlife

CO N FLI C T O F I NTE R E S T
The authors declare that no competing interests exist.

AUTH O R CO NTR I B UTI O N S
MRP and SKG designed the experiment; MRP collected data; MRP, HM, BBS, and SKG analyzed the data and wrote the manuscript.

DATA ACCE SS I B I LIT Y
Data used for ecological analysis are presented in the paper and the supplementary material. Data that identify human participants in the study will not be made publicly available. However, such data can be made available upon request to the authors.