{'fraction_non_alphanumeric': 0.03661342492247157, 'fraction_numerical': 0.02387542098769549, 'mean_word_length': 5.158110882956879, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Soil Physical Quality and Relationship to Changes in Termite Community in Northwestern Colombian Amazon

Conversion from Amazon forest to low-management pasture or agriculture causes not only degradation of aboveground vegetation but also negative changes in soil properties and ecosystem services. This study aimed to evaluate the impact of physical soil degradation on termite community changes in three contrasting land uses (natural regeneration, rubber plantations, and silvopastoral systems). Soil physical quality was assessed through a set of physical variables, such as bulk density, porosity, soil macro-aggregation state, Visual Evaluation of Soil Structure (VESS) and penetration resistance, which were summarized in an overall synthetic indicator of physical quality. Besides, transects of 20 × 2 m were established in each land use; each transect was divided into four sections of 5 m to search and collect termites during 1 hour in each section; likewise, termites were collected from blocks of soil 25 × 25 × 10 cm (length, width, and depth, respectively) adapted from the Tropical Soil Biology and Fertility (TSBF) method. In total, 60 transects were evaluated, 20 in each land use. A total of 41 species were collected across the three land uses evaluated: natural regeneration presented 60% of the collected species (25 species), silvopastoral systems 53% (22 species), and rubber plantations 39% (16 species). Additionally, composition species from the silvopastoral, agroforestry systems, and natural regeneration were different, and a close association between these last land uses was observed. Soil physical characteristics showed significant variations between land uses. The rubber plantations presented lowest values of soil physical quality, while the natural regeneration showed high soil physical quality. These changes affected termite community and lead to changes in its composition with disproportionate loss of some species; however, there are some that can acclimate well to the decline in the soil physical quality.


INTRODUCTION
The Amazon region plays an essential role to the Earth, harboring about 10-15% of land biodiversity (Hubbell et al., 2008;Nobre et al., 2016). It serves as a sink for greenhouse gasses by absorbing up to 5% of annual CO 2 emissions (∼2 billion tons) (Saatchi et al., 2011), and it has been well-recognized as one of the key components of the Earth's climate system (Malhi et al., 2008).
Despite its relevance, this ecosystem is under increasing human pressures related to agriculture and livestock farming, timber extraction, mining, illicit crops, among others, which resulted in unprecedented rates of land cover changes due to forest clearing, degradation, and fragmentation (Sonter et al., 2017;Roque et al., 2018). In Colombia, the Amazon covers 42% of the country's territory (Etter et al., 2006). Despite many efforts to enhance forest landscape conservation and its biodiversity (Furumo, 2020), this region concentrated 70.1% (1,381.8 km 2 year −1 ) of deforestation from a national total in 2018. Deforestation in the department of Caquetá alone accounted for 23.7% (467.6 km 2 year −1 ), becoming an important hotspot of deforestation in Colombian Amazon (SMByC 2020).
Forest clearing not only generates the loss of biological diversity at all levels but also leaves the soil exposed to processes of erosion and degradation. For example, in deforested areas for crops and degraded pastures by domestic livestock, a reduced soil moisture, increase in the bulk density and a decrease in total porosity at soil surface has been observed due to animal grazing, compaction of agricultural machinery, and the impact of the raindrops (Cherubin et al., 2016;de Queiroz et al., 2020).
In tropical rainforests, termites constitute an important part of the soil fauna biomass (Bignell and Eggleton, 2000). Some studies suggest that termites may represent 40-65% of overall soil macrofaunal biomass in some biotopes (Loveridge and Moe, 2004). They contribute to ground decomposition of litter, plant growth, and overall biodiversity (Jouquet et al., 2011;Bottinelli et al., 2015).
Specifically, in tropical and subtropical soils, the role of termites has been widely documented demonstrating its impact on soil structure and properties, such as soil formation and aeration, bioturbation, water content and hydraulic proprieties, soil organic matter, and nutrient cycling (Bignell, 2006;Jouquet et al., 2011;Harit et al., 2017). These alterations are caused by the development of their biogenic structures (nests, foraging tunnels, and formation of galleries) and feeding strategies (ingesting and egesting soil) (Jouquet et al., 2016a).
However, the effect generated by changes in soil characteristics on termite's community distribution has been poorly studied, and it is limited to some efforts aimed to analyze the influence of soil chemical properties on termites (Bourguignon et al., 2015;Jouquet et al., 2015). In this regard, we evaluate the impact of soil physical quality on changes in termite communities and hypothesized that soil physical quality affects different termite communities.

Study Area
The study area was located in the Department of Caquetá southeastern Colombia between 1 • 30 ′ N and 1 • 10 ′ N, and between 75 • 35 ′ W and 76 • 0 ′ W (Figure 1). It has an average yearly temperature of 25 • C, yearly rainfall of 3,600 mm, concentrated between April and November and a light dry season from December to March.
Sampling in this area was carried out during September 2019 in the following land uses: natural forest (NF), rubber plantations (RP), degraded pasture (DP) according to the methodology Corine Land Cover adapted by the Geographic Institute Agustín Codazzi (IGAC) and the Institute of Hydrology, Meteorology and Environmental Studies (IDEAM) for Colombia (Table 1). A brief characterization of the chemical properties of the soil in the evaluated areas is presented in Table 2.

Termite Sampling
Termites were collected from transects of 20 × 2 m (adapted from Swift and Bignell, 2001) located in the central part of each land use type avoiding the edge effects; a minimum distance of 30 m between transects was guaranteed. In the silvopastoral system, the transects were located under the Gliricidia sepium tree line, while in the agroforestry system, they were established under the Hevea brasiliensis and Theobroma cacao intercalated line; a total of 60 transects were evaluated, 20 in each land use.
The transects were divided into four sections of 5 m. In each section, termites were searched for in any specific microhabitat, such as nests, dead wood, trunks, roots, and under rocks, 1 hour per person. When a colony was encountered, collection time could not exceed 3 min to avoid overestimation of the dominant species. Also, termites were collected from extracted blocks of soil by using a 25 × 25 × 10 cm (length, width, and depth, respectively) metallic frame adapted from the International Organization for Standardization (ISO) and the Tropical Soil Biology and Fertility (TSBF) methods (Anderson and Ingram, 1993;ISO, 2011).
Identification was made at genus and species levels, using taxonomic guides and reviews (Constantino, 2002a;Constantino et al., 2006;Constantino and Carvalho, 2012;Rocha et al., 2012;Constantini and Cancello, 2016). When identification was not possible, individuals were separated into morphospecies according to morphological differences. Individuals of the Apicotermitinae subfamily were separated at the species and morphospecies levels by dissecting the enteric valve and comparing its morphology with that described in the literature (Bourguignon et al., 2010(Bourguignon et al., , 2013(Bourguignon et al., , 2016. Termite samples are kept in 2-ml vials containing 80% ethanol, in the Colección Entomológica, Laboratorio de Entomología, Universidad de la Amazonia (LEUA).

Soil Physical Parameters Determination
A set of five physical variables of the soil was evaluated, which was considered key to have a general idea of the physical quality of the soil namely: • Bulk density • Soil moisture • Soil resistance to penetration • Soil macro-aggregation state • Visual evaluation of soil structure Undisturbed soil samples were collected in the center of the 0-10, 10-20, and 20-30 cm layers using a metallic ring (5 × 5 cm 98 cm 2 ) and disturbed samples from the same soil depths. In the laboratory, undisturbed samples were weighed, dried in a forced-air oven at 105 • C for 48 hours, and weighed again. Bulk density (BD, Mg m −3 ) was calculated by dividing the soil dry mass by the volume of the cylinder, whereas soil moisture (%) was determined by the equation: soil moisture = [(dry soil mass/wet soil mass) -1] × 100.
Measurements of Soil Resistance to Penetration (SRP) were performed using a hand penetrometer (Eijkelkamp) around the soil sampling trenches down to 30 cm with angle and surface area of a cone of 60 • and 3.3 cm 2 , respectively. Also, total porosity was calculated by using bulk density and particle density of soil through Equation (1): where TP is the total porosity, Bd and Pd are the bulk density and particle density, respectively. The state of soil macro-aggregation was determined using the methodology described by Velasquez et al. (2007). In each transect, we also collected a small block of soil (10 × 10 × 10 cm) and separated macroaggregates based on the visual separation of macroaggregates (>4 mm) according to biogenic aggregates produced by macroinvertebrates, or root aggregates made of soil stuck to the roots, physical aggregates produced by physical processes, and non-macroaggregated soil. Separation was done by gently breaking the soil apart into its natural constituents.
The Visual Evaluation of Soil Structure (VESS) was performed following the methodology proposed by Ball et al. (2017). In each transect, a soil block of 20 × 10 × 25 cm (width, thickness, and depth, respectively) was extracted, and each block was divided into two layers: topsoil (0-10 cm) and subsoil (10-25 cm). This method involves the removal and gentle breakup of a spadeful of topsoil by hand to reveal the main structural units and any layers of contrasting aggregation. Each layer is compared to the photographs with identified dimensions and descriptions in a colored chart and allocated to one of five soil quality scoring categories (Sq) from Sq 1 = best to Sq 5 = worst topsoil quality (VESS) and Ssq 1 = best to Ssq 5 = worst subsoil quality (SubVESS).
An overall weighted Sq score was calculated for each sample based on the individual score and thickness of each contrasting soil layer, according to Equation (2) VESS Sq score = n i=1 SqiTi TT (2) Land use type Natural regeneration Deforested area that was abandoned ∼15 years ago to recover soil, differing from the forests in species composition (dominated by Cecropia spp.) and in structure (canopy <10 m).

Rubber plantations
Rubber crops planted in late 1990s composed of rubber tree as the main species and some cocoa trees (Theobroma cacao L.) planted since 2015.
Silvopastoral systems An area that combines pastures (Brachiaria brizantha or Brachiaria humidicola) subject to cattle grazing and strips of trees (Gliricidia sepium) isolated and without trampling planted since 2005.
where VESS Sqscore is the overall VESS score of the sample, Sqi and Ti are the score and thickness of each identified soil layer, respectively, and TT is the total thickness of the soil sample.

Termite Community Analysis
The relative abundance of each termite species was expressed as the number of encounters established by recording the presence of each species only once in each plot. Encounters and richness species data per transect were analyzed using a mixed generalized linear model (GLMM) with negative binomial and Poisson distribution for abundance and richness, respectively, stating land use as a fixed effect and farms as a random effect. We use the "lme4" package (Bates et al., 2018). Fisher multiple comparison tests (α = 0.05) were also applied for fixed effect using the "multcomp" package (Hothorn et al., 2008). To compare alphaand beta-diversity across sites, the Shannon exponential (eH ′ ) (Jost, 2006), and Simpson's inversed (1-D) indices, as well as rarefaction and extrapolation curves of richness were calculated using iNext (Hsieh et al., 2020). To test the hypothesis of differences in Shannon exponential and Simpson's inversed between land uses, we performed a generalized linear model (GLM). Non-metric multidimensional scaling (NMDS) based on the Bray-Curtis index of dissimilarity was used to compare composition termite species between land uses; species with a single encounter in all sampling were not considered to avoid the zero-inflated effect. Besides, an analysis of similarities (ANOSIM) was performed with 999 permutations to test statistically whether there were significant differences between land uses. A similarity percentage analysis (SIMPER) was used to identify which species contributed the most to the differences in termite community structure between land use types. These analyses were performed with "vegan" package (Oksanen et al., 2019).

Analysis Between Physical Soil Variables and Land Use Groups
Physical soil variables were analyzed with a generalized linear model (GLM) and a principal component analysis (PCA) associated with a Monte Carlo test from the "ade4" package to assess differences among land use and soil physical characteristics, and their statistical significance (Dray and Dufour, 2007). Besides, physical soil variable dataset was summarized in overall synthetic indicator of physical quality, following the methodology adapted from Velasquez et al. (2007). A principal component analysis (PCA) was used to determine the parameters that best-captured variance within a dataset. We selected those variables with a significant contribution (>50% of the explained variance) in either of the first two main axes. The selected variables were then combined into a single value indicator and scaled to a number ranging from 0.1 to 1.0 using a homothetic transformation.

Analysis Between Physical Soil Variables and Termite Species
To examine the relationship of termite species abundance to physical soil variables, canonical correspondence analysis (CCA) was performed. This analysis was run with two matrices: a dependent matrix containing 24 termite species whose abundance was >1% of the total number of individuals collected, species with abundance lower than 1% were grouped in the category "others" and considered in the analysis. The independent matrix contained the soil physical variables, which we did not keep in the analysis variables that had variance inflation factor (VIF) values higher than 10. The significance of these associations was verified using a permutation test from the "vegan" package (Oksanen et al., 2019); neither were they considered variables that were not significant using permutation tests. All analyses were performed using the R language software, version 3.5.3 (R Core Team, 2018).

Species Richness and Community Structure
A total of 41 species from seven subfamilies and 24 genera were collected across the three land uses evaluated (Table 3) in 135 encounters. Likewise, Nasutitermitinae and Syntermitinae were the most abundant subfamilies with 43 and 24% of the total species collected, respectively. In contrast, Amitermitinae, Coptotermitinae, and Heterotermitinae recorded only one species each in the sampling. The genera with the largest number of species were Nasutitermes, Anoplotermes, and Syntermes containing 10 and 3 morphospecies, respectively, representing 38% of the total specific richness.
There were no statistically significant differences according to GLMM for species richness by transect among land uses; neither were there statistically significant differences according to GLM in the Shannon exponential (H) and Simpson's inverse (1-D) diversity indices per transect between land uses. When analyzing the accumulated species richness in each land use, natural regeneration presented 60% of the collected species (25 species), silvopastoral systems 53% (22 species), and rubber plantations 39% (16 species). Furthermore, sampling coverage curve including interpolation and extrapolation values for each land use averaged 72% of completeness and suggests that the sampling appropriately represent the termite diversity in the study area (Table 4).  Values in bold correspond to families and subfamilies.
The NMDS diagram based on the Bray-Curtis index of dissimilarity showed termite species composition differences (Figure 2). In general, 12 species were found exclusively from natural regeneration, nine from silvopastoral systems, and five in the rubber plantation. Eight termite species were found in all land uses. The ANOSIM demonstrated a significant statistical difference between the land uses evaluated (Global R = 0.29, P = 0.002). The species from the silvopastoral system was different from the species in agroforestry systems and natural regeneration. Moreover, a close association between these last land uses was suggested by the analysis (Figure 2).
The similarity percentage analysis (SIMPER) allowed us to identify which termite species mainly contributed to the dissimilarity in the species composition between land uses. The SIMPER showed that species of the genus Nasutitermes, the xylophages species Heterotermes tenuis, Microcerotermes arboreus, and Nasutitermes guayanae and Aparatermes silvestrii were the main species responsible for the differences between land uses (Table 5). Nasutitermes was the more common species in natural regeneration, with a 12.2% contribution; while Heterotermes tenuis, Microcerotermes arboreus, and Nasutitermes guayanae were the more abundant species in rubber plantation, with a 21.3% contribution, and Aparatermes silvestrii dominant in silvopastoral systems, with a 13.25% contribution.

Physical Soil Variables
Physical soil characteristics showed significant variations between land uses according to GLM. Rubber plantations presented lowest values of the overall synthetic indicator of physical quality (0.41 ± 0.04; Figure 3A) associated with highest values for the parameters of bulk density (1.24 ± 0.06 Mg m −3 ; Figure 3B), soil resistance to penetration (1.17 ± 0.09 Mpa; Figure 3C), and visual evaluation of soil structure (1.49 ± 0.16 Sq score; Figure 3F). Meanwhile, natural regeneration showed high physical soil quality with highest values for soil porosity (57.16 ± 2.69%; Figure 3D) and moisture (44.66 ± 2.28%; Figure 3E).
The variance of physical soil properties was explained by 21.24% by land use according to the PCA. Axes 1 and 2 of the PCA represent 36.2 and 19.2%, respectively. Axis 1 separated on one side the sites with high porosity, macro-aggregation status, and soil moisture ( Figure 4A) to which secondary vegetation was associated, at the other end of Axis 1 are the sites with high values for the parameters of bulk density and soil resistance to penetration like rubber plantations ( Figure 4B). The formation of Axis 2 was mainly contributed on one side by sites with high remaining material from aggregate morphology and on the other side by high Sq score visual evaluation of the soil structure.
Soil macro-aggregation on average, 21.5% of the soil volume was not macro-aggregated, biogenic aggregates contributing 21.8% of the total soil mass and physical aggregates comprising 20.9% (Figure 5). Important differences between land uses were suggested in soil macro-aggregation, with high fraction from remaining material representing 33.4% of soil under rubber plantations, while natural regeneration possessed much larger percentages of biogenic aggregates (28.1%).

Effect of Physical Soil Variables on Termite Species Abundance
The CCA explained 16.1% of the variation in termite species (P < 0.001) and included variables related to soil physical variables (bulk density, soil resistance to penetration, total porosity, visual evaluation of soil structure, and soil moisture). The first two axes of the biplot accounted for 48.9 and 45.0%, respectively, of the variance explained by the CCA (Table 6). Some termite species were associated with the physical soil variables, for example, Embiratermes neotenicus associated with high soil moisture, Curvitermes odontognathus and Cyrilliotermes angulariceps with high total porosity, a species group made up of the species Cylindrotermes capixaba, Microcerotermes arboreus, and Termes hispaniolae was related to high values of visual evaluation of soil structure, finally soil resistance to penetration was correlated positively with some species of the genus Nasutitermes and negatively with total species richness (Figure 6).

Termite Community
Although the number of species per transect did not present significant differences, the total number of species did show important variations, the area of natural regeneration to be on the way of recovery of termite diversity with 25 species collected  in total, if it is compared with previous studies carried out in the Amazon that report 52 species for forest (Ackerman et al., 2009). The species richness found in the rubber plantation in our study (16 species) was greater compared to others' studies in Colombia for commercial rubber plantation (10 species; Pinzon et al., 2012). Despite the present study was not evaluated, in the area of natural forest, the total diversity observed was similar to studies carried out by Pinzon et al. (2017) who reported the presence of 40 termite species in riparian forests of the eastern plains of Colombia. However, these data are not completely comparable because of the differences between land uses and the study areas; it can be used as a reference.
Several authors have reported a negative effect on termite assemblages according to the level of habitat disturbance and FIGURE 4 | Projection of physical soil characteristics and the 60 sampling points according to land uses (A), and in the plane formed by Axes 1 and 2 of PCA (B) in Colombian Amazon region. Ellipses show sampling points grouped from a given land use centered on their respective bari-centers (P < 0.01: Monte Carlo permutation test, Explained variance: 21.24%). BA, biogenic aggregates; BD, bulk density; NA, soil not macro-aggregate; PA, physical aggregates; RA, root aggregates; RM, remaining material; SRP, soil resistance to penetration; TP, total porosity; VESS, visual evaluation of soil structure. trampling, with a reduction in termite diversity and abundance (Barros et al., 2004;Vasconcellos et al., 2010). Nevertheless, the inclusion of trees in pastures and the isolation of these to prevent livestock from entering have benefited the diversity of termites in our study site (22 species) exceeding the number of species reported for degraded pastures (Carrijo et al., 2009;Duran-Bautista et al., 2020). The changes in species composition observed in this study confirm the results of previous work that suggests that simplified habitats are associated with a decrease in termite richness (Vasconcellos et al., 2008;Luke et al., 2014). Besides, the exclusive presence and high dominance of species, such as Microcerotermes arboreus in the rubber plantations, suggest the creation of favorable environments that could turn them into pests (Constantino, 2002b). Likewise, in Colombia, the Hetereotermes spp. have been reported in high densities in rubber crops of Puerto López-Meta (Pinzon et al., 2012) and as a pest in rubber plantations in the Amazon region (Sterling et al., 2011).
Previous work has shown that Nasutitermes spp are frequent in areas with regeneration process for 5 years (Viana-Junior et al., 2014) as well as in reforestation sites with native species (de Paula et al., 2016).

Land Use Effect on Physical Soil Variables
In general, all land uses presented a good soil physical quality, and the natural regeneration presented the best values. In this sense, it is proposed that changes in forest land use to more intense land uses, such as agroforestry or silvicultural systems generate changes due to disturbance that alters the chemical and physical composition of the soil (Tellen and Yerima, 2018). The time of abandonment and non-intervention in the areas of natural regeneration seems to have generated a positive effect that lead to soil physical characteristic protection (Lisboa et al., 2009), reflected in high values of total porosity and soil moisture, low resistance to penetration, and consequently better overall quality of physical soil characteristics.
However, important differences were observed in rubber plantations with higher values of apparent density and resistance to penetration even surpassing silvopastoral systems. This could be explained by two mechanisms: (1) the silvopastoral systems evaluated were recently established, and (2) the samples were taken below the tree line where the cattle did not have access. de Souza Braz et al. (2013) propose that in areas where grazing pressure is low, it avoids the negative effects on physical soil characteristics generated by cattle trampling, such as increasing bulk density and soil resistance to penetration, reducing soil aeration, aggregate stability, and water infiltration that have been widely documented (Arevalo et al., 1998;Sharrow, 2007;Drewry et al., 2008).
Concerning the above, a study found that rubber agroforestry systems with low diversification, present higher values of apparent density and resistance to penetration than more diversified agroforestry systems due to high abundance of roots mixed with a soft layer composed primarily of porous and rounded aggregates present in these systems (Cherubin et al., 2018). Likewise, an assessment of the effect of different land uses on soil quality suggests that rubber trees are the plant cover causing the smallest change in terms of soil quality (Moreira et al., 2011).

Physical Soil Variables and Relationship With Termite Species
Soil compaction indicated by resistance to penetration had a negative relationship with total termite species richness. In this regard, Jones et al. (2003) suggested that compaction of soil surface and an increase in soil bulk density might have had a negative impact on soil-dwelling termites. This may be due to creation of tunnel and foraging galleries during termite's feeding and foraging activities by moving soil particles of different sizes (Jouquet et al., 2011(Jouquet et al., , 2016b. These activities are difficult by increased soil compaction. Tucker et al. (2004) found that construction speed of tunnel network varied with soil compaction, being much slower in the more compacted soil (1.35 g/cm 3 ). Although this study was conducted under laboratory conditions using a different substrate and moisture content, and it cannot be directly compared with the present study, however, clearly the same relationship was observed.
Nevertheless, some species of the genus Nasutitermes manage to support well in conditions of high resistance to penetration, and this would be related to the ability of these species to adapt to disturbed environments. For example, they can be pests in sugarcane crops (Constantino, 2002b;Miranda et al., 2004) where the physical quality of the soil is low by compaction due to intensive mechanization (Cherubin et al., 2016). The relationship of E. tenicus with high-humidity soil could be explained by two mechanisms: the first is the nesting habits developed for this specie as it builds an epigeal earthen nest of variable shape, sometimes at the bases of trees (Constantino, 1992a,b), and the second is one of the more frequently found termite species in primary and secondary forests (Bandeira et al., 2003;Davies et al., 2003;Souza et al., 2012) where soil moisture is higher.
Cyrilliotermes angulariceps was positively related to total porosity; the species of this genus are soil feeders (Constantino and Carvalho, 2012), and specifically, C. angulariceps feeds on soil particles richer in organic matter (Bourguignon et al., 2009) where soil porosity is generally higher. Although little information is available on the influence of termite foraging activity on soil porosity (Jouquet et al., 2016a), some authors have demonstrated that termite mediated increases in soil porosity (Lavelle, 1997). The relation of Cylindrotermes capixaba, Microcerotermes arboreus, and Termes hispaniolae with high values of soil structure can be explained because these species were presented in high density in the rubber plantations where this parameter presented the highest values. These results coincide with those from Thoumazeau et al. (2019), which reported a VESS gradient from intensive cash crops, the most disturbed system, through different ages of rubber planting to forest as the least disturbed system.
In general, changes in physical soil quality presented in different land uses affect termite communities and lead to changes in its composition with disproportionate loss of some species; however, there are some who can adapt well to the decline in physical soil quality, especially in rubber plantations and silvopastoral systems. Among the parameters evaluated, the resistance to penetration was the one that generated the greatest negative impact, and soil moisture was positively related to the presence of some species.

DATA AVAILABILITY STATEMENT
The datasets presented in this article are not readily available because they were financed by the University of the Amazon.
Requests to access the datasets should be directed to Ervin Humprey Duran-Bautista, ervinduranb@gmail.com.

AUTHOR CONTRIBUTIONS
ED-B: conceiving and designing analyses, data collection, data analysis, and wrote paper. YM, JG, TO, and MB: data collection, data analysis, and wrote paper. All authors contributed to the article and approved the submitted version.