Introduction

Woody plants in African savannas provide beneficial services for human livelihoods including timber and fuelwood, food sources for domestic and wild browsers and carbon storage (Makhado et al., 2014). While woody plants are integral components of savanna ecosystems, tree and shrub cover has increased over the past century in savannas worldwide, a process known as bush encroachment (Venter et al., 2018). Bush encroachment often results in negative economic and ecological shifts including decreased grass biomass used for livestock production (Archer and Predick 2014), decreased biodiversity (Ratajczak et al., 2012), and depleted soil moisture via increased rainfall interception and transpiration (Honda & Durigan, 2016). In southern African savannas, bush encroachment is associated with a shift from tall trees to short-statured shrubs (Zhou et al., 2021), which can make it difficult to see the charismatic animals that attract tourists, possibly resulting in a loss of tourism revenue for protected areas (Gray & Bond, 2013;Luvuno et al., 2022).

Frequent or high-intensity fire is often proposed as the most effective way to prevent largescale bush encroachment (Scholtz et al., 2022), but the benefits of high-intensity fire to reduce woody cover are short-lived (Strydom et al., 2022). Additionally, in semi-arid and bush encroached savannas, grass biomass is often too low to carry regular fires that are intense enough to limit woody cover (Smit et al., 2013;Scholtz et al., 2022). More intensive management approaches are needed in these systems to control woody plants, open the canopy, and increase grass biomass. Bush clearing has become a common management technique to mitigate the negative effects of bush encroachment (Smit, 2005;Eldridge & Ding, 2021), particularly in small game reserves where fire is a risk to infrastructure (Smit et al., 2022). This strategy typically involves the removal of small-statured trees and shrubs while maintaining large trees to preserve habitat heterogeneity, shade for livestock and wildlife, and the gradients of nutrient and water availability associated with large trees (Schmitt et al., 2022). Bush clearing is labor-intensive and expensive (Luvuno et al., 2022), but may be a viable and necessary alternative in areas where fire is infrequent and herbicide application is undesirable.

Long-term reductions in tree and shrub cover require either the use of herbicide or repeated clearing efforts, as most savanna woody species are well-adapted to disturbance and resprout following the loss of their above-ground parts (Bond & Midgley, 2001). These species use nonstructural carbohydrates (NSC) stored in their trunks and belowground organs to regrow aboveground tissues following disturbance (Wigley et al., 2019). Single or infrequent disturbances only have short-term effects on tree size and cover because these species typically have sufficient NSC stores to recover. In addition, resprouting often results in the transformation of tall, single-stemmed trees into short, multi-stemmed shrubs which can make game viewing more difficult and contradict the initial goals of bush clearing efforts (Fisher et al., 2014).

Repeated cutting has the potential to deplete NSC reserves in woody species to the point of mortality -especially in conjunction with other pressures like browsing or fire -leading to longterm reductions in woody plant cover (Smit, 2004).

Our focus was on Colophospermum mopane (Kirk ex Benth.) J. Léonard (hereafter 'mopane') a dominant tall tree or short multi-stemmed shrub (< 4 m) that occurs over ~ 35% of Southern African savannas (Stevens, 2021). The short-statured shrubs have increased in abundance over the past several decades (Zhou et al. 2021) resulting in encroachment (O'Connor et al., 2014;Stevens, 2021). Given the widespread dominance of mopane across Southern Africa, it is necessary to assess the potential costs and benefits of clearing this species, particularly in protected areas as most research on the effects of bush clearing have been conducted in rangelands (Smit 2004;Fisher et al., 2014).

Here, we used repeated bush clearing to ascertain if this treatment resulted in wide-spread shrub mortality, and describe any long-term changes in ecosystem structure including changes in grass biomass, herbivore presence, and woody plant physiology. In this experiment, mopane shrubs < 4 m tall were cut repeatedly from 2015-2022 to explore if repeated bush clearing is a viable alternative to herbicide and fire. We assessed the effects of cutting on grass biomass and monitored herbivore presence via dung counts. We predicted that repeated clearing would increase grass biomass, due to reduced competition between mopane and grasses, and increase grazer abundances, due to increased forage and a preference for more open areas by certain herbivore species (Burkepile et al., 2013). We also examined the effects of cutting on mopane belowground NSC storage and functional rooting depth, to assess the mechanisms by which repeated cutting eventually leads to mortality and long-term alteration of woody plant abundance. We predicted that repeated bush clearing would deplete belowground NSC storage, reducing the energy available for resprouting and shift mopane water uptake to shallower soils due to potentially reduced carbon allocation in deeper roots (Landhäusser & Lieffers, 2011).

Methods

Study site -Mthimkhulu Game Reserve is a rural 7500-hectare community-owned reserve that shares an open border with Kruger National Park, South Africa (23°31'46' S; 31°06'12' E) and is managed by the Mthimkhulu Tribal Authority. Since 2016, the South African Environmental Observation Network has been working with this authority on long-term, socio-ecological research that can both improve ecosystem services and create employment for rural people living adjacent to the reserve. The tribal authority identified eco-tourism as the most desirable land use for the reserve, and bush clearing was suggested to enhance tourism, as has typically been done in privately-owned game reserves in the region.

The site receives 467 mm mean annual rainfall, most of which falls between November and April (Mahlangeni Ranger station, 14 km south of the site, 1968 -2010). Summers are hot and humid, with maximum temperatures typically above 30 °C, while minimum temperatures in winter rarely fall below 10 °C. The region has granite-derived, nutrient poor soils but the experimental plots were located on an old alluvial terrace of the Klein Letaba River with soils at least 1 m deep. Soil texture is approximately 62% sand, 20% clay, and 17% silt (Table S1). The site falls within the Lowveld Mopaneveld vegetation type of Mucina and Rutherford (2006), a semi-arid savanna characterized by a dense cover of mopane shrubs, sparsely scattered trees, and low grass cover. The plots contained dense monotypic mopane stands with occasional Combretum imberbe Wawra and Vachellia tortilis (Forssk.) Gallaso & Banfi trees and Grewia bicolor Juss. shrubs. The grass layer was sparse and consisted mostly of annual Aristida species, with scattered tufts of perennial grasses Urochloa mosambicensis (Hack.) Dandy and Panicum maximum Jacq. Most of the ground cover consisted of bare soil and mopane leaf litter. Fire is infrequent in the region (~10 years; Smit et al., 2013), but the site burned in 2014, the year before the study began, resulting in top-kill of most mopane shrubs <2 m (personal observation).

Experimental design and community sampling -In 2015, we established six 60 x 60 m plots divided into three blocks (two plots per block). The blocks were spread equidistant across an alluvial terrace, with specific locations selected to ensure that tree and shrub cover was similar within each block. One plot within each block was randomly assigned to a repeated mopane clearing treatment and the other plot assigned as a control (n = 3 cleared and 3 control plots).

From 2015-2022, mopane shrubs < 4 m tall were cut at the base 2-3 times during each growing season with a panga and cut material was removed from the plots (Fig. 1). The timing of each clearing was determined by the height of resprouting stems, with clearing taking place shortly after the start of each growing season (October -December) and then again once resprouting stems reached 30 -50 cm tall. For most response variables, sampling began in 2017 at the end of the growing season. In every plot, volumetric soil water content (VWC) was measured every hour (Oct. 18, 2017 -June 13, 2022) using CS655 time-domain reflectometry probes (Campbell Scientific Inc., Logan, UT, USA) at 10, 30 and 80 cm depths. To assess the effects of mopane clearing on the plant and animal community, we randomly established two 50 m transects within each plot. Transects were nonoverlapping and at least 20 m apart and 10 m from the nearest plot boundary. Grass and shrub cover were estimated using the point-step method (Evans & Love, 1957). At each 1 m point along each transect we measured the distance to the nearest perennial grass tuft, the width of the tuft, and the presence of mopane cover directly above each 1 m point.

Mopane individuals were categorized into four height classes: < 0.5 m, 0.5-2 m, 2-5 m, and > 5 m. The number of mopane in each height class was summed across the two transects in each plot. Frequency of each height class was calculated as the percentage of points with tree or shrub cover directly overhead. At the end of each growing season, herbaceous biomass (g m -2 ) was clipped in four -six 0.5 x 0.5 m quadrats randomly located per plot. Biomass samples were dried at 80 °C for two days, sorted into current and past season material, and weighed. We used each transect line as a 4 m wide belt transect to estimate animal visitation by counting species dung piles approximately every 1-2 months from April 2016 -August 2022.

Mopane physiology-At the end of each growing season, mopane mortality was surveyed in each cleared plot by walking from a random cut mopane to its nearest dead or alive neighbor until 100-250 individuals were recorded. Individuals were considered dead if no green resprouting tissue was present. As the stumps of dead individuals persisted for many years, this provided a viable means of estimating mortality over time. We did not measure mortality in the control plots, but background mortality rates in mopane-dominated savannas are typically below 5% (Swemmer, 2020).

We used established protocols to measure functional rooting depth of mopane and adjacent grasses using the stable isotopic signature of xylem-and soil-water (Dawson et al., 2002). In March of 2015, 2017, and 2018, we collected non-photosynthetic tissue from mopane (stem tissue) and grasses (crown tissue) within each plot to collect xylem water. We collected soil samples from 10, 30, and 50 cm depths from soil pits. All samples were placed in exetainer vials on ice, and frozen until processing. We used cryogenic vacuum distillation to extract water for isotopic analysis. Samples were analyzed on the Picarro L1102-I CRDS analyzer (Picarro, Inc., Santa Clara, CA) at Kansas State University. Samples were referenced to V-SMOW and converted to delta notation using:

Where R is the ratio of the heavy to light isotope. ChemCorrect software was used to flag samples that had organic contamination, which were removed from any subsequent analyses.

We measured NSC concentrations of mopane belowground stems to assess changes in carbon storage. Sections of mopane stems/boles were collected 10-20 cm below the surface concurrent with collection for isotope samples. Samples were washed to remove soil, microwaved for 90 sec to stop enzymatic activity, dried for 72 h, and ground to a fine powder.

NSC extraction and analysis followed the procedure outlined in O' Connor et al. (2020).

Statistical analysis -We used R V4.2.1 (R Core Team, 2022) for all statistical analyses.

Repeated measures ANOVAs using the lme4 and car packages (Bates et al., 2015;Fox & Weisburg, 2019) were used to test the effects of mopane clearing on soil VWC, mopane and grass responses, and mopane NSC concentrations. We used Tukey's HSD for pairwise comparisons using the package emmeans (Lenth, 2022). For soil VWC, distance to perennial grass tuft, width of grass tuft, and mopane NSC analyses we included an interaction between treatment and year as fixed effects and block as a random effect. For soil VWC, we only included data during the growing season (November -April). We ran the same model for grass biomass, but included quadrat nested within block as a random effect to account for variability due to multiple quadrats taken within each plot (n = 4-6 quadrats per plot). For the frequency of mopane size classes, we included an interaction between treatment and size class as fixed effects, with block nested within year as random effects. To meet the assumptions of normality and homogeneity of variance, we transformed grass biomass (square root) and mopane frequency, distance to the nearest perennial grass tuft, and the width of the grass tuft (log+1).

Dung counts were analyzed following Voysey et al. (2021). Since the length of transects used for dung counts varied from 50 -100 m, we standardized dung counts by transect length and summed the number of dung piles per species across the two transects within each plot. We then divided the number of dung piles by the number of days between sampling dates to estimate dung deposition per day. To assess if this dung metric significantly differed for each species between cleared and control plots, we used a generalized mixed effects model with a Tweedie distribution using the glmmTMB function with a log link function and date nested within block as a random effect (Brooks et al., 2017). The Tweedie distribution was used to account for zeroinflated data common in herbivore dung surveys (Voysey et al., 2021). We calculated herbivore preference for cleared vs. control plots by dividing the number of dung piles in cleared plots by the total number of dung piles in cleared and control plots for each sampling date, where 1 is complete preference for cleared plots and 0 is complete preference for control plots (Voysey et al., 2021;Donaldson et al., 2018).

To assess differences in functional rooting depth among grasses and mopane in the cleared and control treatments, we analyzed each year (2015, 2017, and 2018) separately. Since 𝛿𝛿D and 𝛿𝛿 18 O were collinear and varied similarly with soil depth (Figs S1-S2), we collapsed 𝛿𝛿D and 𝛿𝛿 18 O into a single axis using PCA to analyze both water isotopes in a single analysis (sensu Holdo et al., 2018;Case et al., 2020). The PCA approach is an alternative to isotope mixing models when plant signatures fall outside of the range of sampled soil signatures, indicating plants are using deeper soil water than sampled. We then used PC1 as the response variable in a linear mixed effects model with an interaction between vegetation type (mopane or grass) and treatment and included block as a random effect.

Results

Soil moisture -VWC was higher in cleared than control plots at all depths (P < 0.001; Fig. 2a; Table S2), but these differences were first evident at 10 cm soil depth. Across all growing seasons, mean daily VWC was 39% higher at 10 cm, 10% at 30 cm, and 18% at 80 cm.

From 2020-2021, when rainfall was above-average (Fig. 2b), the difference at 10 cm was nearly 60%.

Vegetation responses -Mopane mortality was low in the first few years of clearing, then increased steadily to 80% cumulative mortality after 7 years (Table 1). The frequency of mopane taller than 0.5 m was significantly lower in cleared than control plots (Fig. 3; Table S3). In cleared plots, most mopane were less than 0.5 m tall, while most mopane in control plots were between 2-5 m tall with few small trees.

Grass biomass was higher in cleared than control plots in 2018 and 2020-2022 (Fig. 4a; Table S4). For these years, average grass biomass was 1.7, 2.1, 2.2, and 1.4 times higher in the cleared than control plots, respectively. The mean distance to the nearest perennial grass tuft was shorter in the cleared treatment in 2021, suggesting higher average grass density in cleared plots (Fig. 4b). The widths of perennial grass tufts did not differ among treatments, except in 2017, where tufts in the cleared plots were larger than those in the control plots (Fig. 4c).

Animal dung counts -Between April 2016 and August 2022, dung from 11 herbivore species were recorded (Table S5). Dung of impala (Aepyceros melampus), buffalo (Syncerus caffer), elephant (Loxodonta africana), and giraffe (Giraffa camelopardalis) were the most common and present in all years, while many herbivores had low dung counts and were not present in all years (Fig. S3). Total dung counts across years were higher in cleared than control plots, mostly driven by buffalo and impala (Tables S5, S6). These species had significantly higher dung deposition per day in the cleared plots than control plots (Fig. 5; Table S6) and preferred cleared plots to control plots (Fig. 6).

Mopane physiology -Glucose and sucrose concentrations were unaffected by cutting in most years except for 2017 (Fig 7a, b). In contrast, stem starch concentrations were significantly lower in cut than control mopane in 2016-2018 (Fig. 7c; Table S7). This difference increased with repeated cutting, with average starch concentrations 1.8, 2.6, and 7.8 times greater in control than cut shrubs in 2016, 2017 and 2018, respectively.

Soil water isotopic composition represented by PC1 declined with depth in all years

where shallow soil layers had isotopically enriched signatures (Fig. S4). Mopane used deeper soil water than grasses in cleared and control plots (Fig. 8; Table S8). In 2017 and 2018, cleared mopane shrubs tended to use shallower soil water than control shrubs, but these differences were only significant in 2017 (Fig. 8).

Discussion

Repeated clearing of a dominant woody species initially appeared to have only minor effects on mortality (5% after 2-3 years of clearing). However, repeated clearing eventually resulted in high shrub mortality (cumulatively 80% after 6-7 years), with increased soil water availability, greater grass biomass, and altered herbivore activity. This experiment demonstrates that repeated cutting can be an effective management tool for opening the canopy after bush encroachment, at least in mopane-dominated, semi-arid savannas. While some ecological responses occurred rapidly (after only a few cuts), multiple cuttings over many years were required to induce significant mortality and long-term reductions in woody cover.

Grass responses -Woody plants in semi-arid savannas have been hypothesized to have facilitative effects on grass growth and biodiversity through amelioration of heat and water stress and nutrient enrichment under canopies (Dohn et al., 2013;Moustakas et al., 2013). However, woody plants in high densities can reduce grass biomass and diversity through shading (Pilon et al., 2020). Our results indicate that mopane trees have strong competitive effect on the grass layer, and that reducing water in the upper soil layers is likely to be a mechanism of competition.

Higher woody cover can lead to decreased moisture at surface and deep soil layers through increased transpiration, canopy interception of rainfall, and reduced infiltration (Smit & Rethman, 2000;Aldworth et al., 2022). There appeared to be a lagged response of 5 years after initial clearing before grass biomass consistently increased in the cleared plots, although this was confounded with differences in annual rainfall (the first three years of the study period had below-average rainfall, and the last three years above-average). Changes in perennial grass tuft size and 'point-to-tuft' distance over time suggest that while perennial tufts in the cleared treatment were able to increase in size after the initial cuts, increased density occurred after two years of cutting. This increase in density presumably created substantial increases in biomass over the entire sward. The large interannual fluctuations in 'point-to-tuft' distances suggest that the dominant perennial species (i.e., Urochloa mossambicensis) had population fluctuations in response to changes in annual rainfall, but with greater recruitment in the cleared plots. This is consistent with the results of Wilcox et al (2020) who found a rapid recolonization by this species following a severe drought at a nearby site.

Herbivore Responses -The opening of the canopy that facilitated increased grass biomass also created areas of preferred habitat for impala and buffalo. We lack the required data to mechanistically explain habitat selection, but results could be due to a combination of increased forage quantity, grass species palatability, and predator visibility (Burkepile et al., 2013). Our data suggest the 60 x 60 m plot size was meaningful for impala and buffalo abundance, but other megaherbivores, such as elephants, are unlikely to show preferences at this scale due to their movement across large areas. For other species, low dung counts also limited our ability to draw general conclusions on herbivore feeding guilds (Table S5). Further research is needed to test whether these patterns would hold if clearing were implemented over larger scales. At the landscape scale, it is likely that the grazer density would either be unaffected or increase, given the greater supply of forage and shade availability provided by sporadic woody cover. Browser density could decline if vast areas of woody cover were cleared. Ideally, reserve management would aim to create a mosaic of cleared and uncleared landscapes, similar to Schmitt et al. (2022), who found herbivore species richness and abundance were highest in areas of intermediate woody cover, and lowest in areas of low habitat heterogeneity (similar to the monodominant mopane stands in the control plots of this study).

Mopane physiology -Despite repeated loss of all above-ground parts, most mopane individuals persisted for many years. This demonstrates the remarkable resprouting ability of mopane and highlights the difficulty of creating long-term reductions in shrub density. However, reduced starch reserves weaken the ability of mopane shrubs to recover following disturbance and the combination of multiple disturbances (e.g, cutting with drought, fire, or herbivory) may increase shrub mortality and long-term, large-scale reductions of woody cover (O'Connor et al 2020; Staver et al., 2009). Repeated defoliation can reduce carbon investment in roots (Landhäusser & Lieffers, 2011), which may have driven the shift in water-uptake to shallower soil layers in cut mopane shrubs. However, cleared mopane still used deeper soil water than grasses, suggesting that deeper soil water is critical for their survival and hydrological niche partitioning exists even when shrubs have been depleted of belowground reserves. Mopane invests most of its fine roots in in the top 40 cm of soil but coarse roots between 40-60 cm of soil (Smit & Rethman, 1998). This shallow and sprawling root system suggests differences in water uptake depth between trees and grasses may be small, particularly when water is available in the surface layers (Kulmatiski et al., 2020).

Large-scale and long-term changes -Transition of a savanna ecosystem from a bush encroached to a stable open state requires internal reinforcing feedbacks after ceasing tree removal treatments (Ratajczak et al., 2018). Scholes (1990) found that mopane stands cut and treated with herbicide returned to their pre-cleared density within 14 years, suggesting that onetime bush clearing is unlikely to reduce mopane cover over decadal time scales, and clearing is required about every 10 years. However, in that study, the application of herbicide was estimated to cause only 40% mortality. We found 80% mortality after 7 years of repeated cutting and return to original densities may therefore take much longer. While the method used for this study is labour-intensive and requires greater short-term investment, it may prove to be more costeffective in the long term. Rates of seedling recruitment will play a key role in determining effective control, as seed sources from uncleared areas could facilitate large recruitment events in years with ideal environmental conditions. Clearing larger areas, with lower perimeter:area ratios could be beneficial in reducing seed establishment from surrounding areas. Additionally, recruitment events may be limited by increases in grass biomass that limits woody plant establishment and growth rates and supports more intense fire, particularly in wet years when grass biomass is high (February et al., 2013;Riginos, 2009).

The economic value of bush clearing encroached savannas in South Africa has been estimated to be US $2.1 billion (Stafford et al., 2017). This estimate included the potential benefits of restoring soil water recharge and grazing capacity and using harvested wood for electricity, fuel, and wood composite products. This presents potential industrial economic benefits of clearing bush encroached areas on a large scale. Additionally, bush clearing may benefit local individuals at smaller scales by using targeted woody biomass for firewood and fencing (Makhado et al., 2014). However, we recognize bush clearing has variable effects dependent on climate and the woody species to be removed (Eldridge & Ding, 2021), and proper ecosystem management will vary depending on the goals of local communities, ecologists, and economists. The method used in this study was notably more expensive (~US $22,000 ha -1 for seven years of repeat cutting or ~US $3,100 ha -1 yr -1 ) than the widely used method of a one-off cut followed by herbicide application (~US $4,300 ha -1 for one year). However, the use of herbicides may be undesirable due to the toxicity of the chemicals to humans and wildlife and the potential effects on non-target plant species. Furthermore, our method created regular employment in an area where poverty and unemployment are high. While repeated clearing provides a limited amount of temporary and unskilled labour, similar interventions at large scales and the resulting transformation of dense savanna to open savanna could increase ecotourism potential. Increased ecotourism is a primary goal of the game reserve used in this study and could create long-term, skilled, and higher paid employment opportunities in areas where few economic opportunities exist. The social benefits may outweigh financial costs, particularly if donor funding can be used to support management objectives. Finally, more intensive clearing methods may be necessary in the future, if the drivers of encroachment, such as elevated atmospheric CO2, continue to increase. Targeted and frequent management practices are increasingly necessary in a changing climate and fragmented landscapes where the return of disturbance (fire or browsers) to the system is insufficient to limit bush encroachment (Case & Staver, 2017;Collins et al., 2021).

In summary, our results highlight that repeated, targeted bush clearing can increase soil moisture, grass biomass, and ungulate habitat in semi-arid savannas. Crucially, repeated clearing was necessary to reduce starch reserves to levels that prevent resprouting and induce tree mortality. While this approach was more expensive than conventional methods, it had positive ecological effects, created social benefits in the form of temporary employment, and may produce longer-lasting reductions in shrub densities. We recommend decadal-scale management plans for bush-clearing projects that include repeated clearing to increase the return on investment and successfully mitigate the negative impacts of bush encroachment in the long term. More research is required to determine the exact thresholds of shrub mortality and understand how less intensive bush clearing in combination with other disturbances, such as fire, may achieve similar management goals dependent on the site.

Statement on Inclusion

Our study took place on Mthimkhulu Game Reserve, a community-owned game reserve in Limpopo, South Africa. Residents from Phalaubeni village, adjacent to the reserve, were employed as labourers for bush-clearing and assisted with data collection. Community members were involved in the planning and execution of this experiment. monodominant stands. (b) Plots cleared of all mopane trees and shrubs < 4 m tall. (c) resprouting mopane shrub. Photos taken in March 2020. and 80 cm soil depths for November 2017 -July 2022. Growing season spans November -April. (b) Total monthly precipitation (mm) from November 2017 -July 2022. Inset shows total precipitation (mm) for each hydrological year (July 1 -June 30). 2022). Bush clearing began in 2015. The height of every mopane tree at each meter along two 50 m long transects was measured and categorized into one of four size classes: < 0.5 m, 0.5-2 m, 2-5 m, or > 5 m. The number of mopane in each height class was summed across transects and frequency was calculated as the percentage of points with mopane cover directly overhead. ** p < 0.01, *** p < 0.001. grass tuft (cm), and (c) width of nearest perennial grass tuft (cm) in cleared and control plots. * P < 0.05, ** P < 0.01. Grass data from 2020 and 2022 are missing due to sampling constraints. Species are ordered from the highest total dung counts recorded to the lowest. Standardized dung counts were averaged for each species in each plot across all sampling dates (n = 3). Note panels (b-d) are an order of magnitude higher than panels (e-i). Panel (a) has a unique scale. *** P < 0.001. calculated from the total standardized dung counts for each species in each plot across all sampling dates. Preference ratio ranges between 0 and 1, where 1 is complete preference for cleared plots and 0 is complete preference for control plots. Grey points are the preference ratio for each plot (n = 3). mopane trees in cleared and control plots. ** P < 0.01, *** P < 0.001. soil water. Asterisk represents significant differences between cleared and control mopane trees in 2017 (P < 0.05).