Restoration thinning permits stems to capitalize on high-rainfall years in a regenerating endangered forest ecosystem

1. Passively regenerating native vegetation presents a cost-effective opportunity to sequester carbon and reinstate habitat in heavily cleared agricultural landscapes. species interactive effects of neighbourhood and moisture availability on the and of A. harpophylla Our clear thinning A. harpophylla interspecific and cristata is relatively weak. As such, thinning of dense harpophylla be with seeding or planting of co-occurring tree species with complementary niches to further accelerate forest recovery in this extensive regrowth ecosystem.

to the scale of historical clearing. By contrast, passively regenerating native vegetation presents cost-effective opportunities to sequester carbon and reinstate habitat in modified landscapes (Dwyer, Fensham, Butler, & Buckley, 2009, Fensham & Guymer, 2009, Geddes, Lunt, Smallbone, & Morgan, 2011. However, passively regenerating vegetation on cleared land may not follow the same successional pathways observed under natural disturbance regimes. Prolific recolonization by one or a few woody species can occur resulting in dense, low-diversity stands that are slowgrowing and slow to self-thin (Jones et al., 2015, Vesk & Dorrough, 2006. While these 'dense regrowth' or 'dog hair stands' (Baskin, 1999, McHenry, Wilson, Lemon, Donnelly, & Growns, 2006) may eventually self-thin, selective thinning of over-dominant species (restoration thinning) has been proposed to reduce intraspecific competition, and accelerate stand development towards a mature state (Stone, Kolb, & Covington, 1999, Swinfield, Afriandi, Antoni, & Harrison, 2016. Thinning of over-dominant species may also accelerate the growth of co-occurring species that were formerly common, as has been widely documented in the forestry literature (Oliver & Larson, 1996), but the magnitude of this effect depends on how the co-occurring species compete for limited resources (Adler et al., 2018).
Restoration thinning has been applied to a wide variety of regenerating ecosystems to meet a range of objectives. For example, restoration thinning in dense Eucalyptus regrowth in Australia has been shown to increase native understory plant cover and species richness (Brown, Murphy, Fanson, & Tolsma, 2019, Jones et al., 2015 and native reptile richness and abundance (Craig et al., 2010). Restoration thinning may also deliver the co-benefit of accelerating stand-level carbon sequestration rates (Dwyer, Fensham, & Buckley, 2010a, Vargas, Allen, & Allen, 2009).
One of the most extensive regrowth ecosystems in Australia is brigalow, which, in its mature state, is characterized by the presence of the brigalow tree (Acacia harpophylla). Depending on the location and soil type, brigalow forms associations with belah (Casuarina cristata), yellowwood (Terminalia oblongata), Eucalyptus species and species from semi-evergreen vine thickets (e.g. Brachychiton rupestris). Prior to European settlement, brigalow forests and woodlands occupied more than 7 million ha of inland Queensland and New South Wales (in the region now known as the Brigalow Belt). In particular, brigalow was dominant on fertile clay soils between the 500 and 750 mm annual rainfall isohyets (Isbell 1962). Less than 10% remains in a mature state due to clearing for pastoral production and dryland agriculture (Accad, 2001).
However, A. harpophylla can prolifically resprout on cleared land via root suckers, and more than 700,000 ha of regrowth persists throughout the Bioregion (Lucas et al., 2014). While dense regrowth stands eventually self-thin over many decades (Johnson, McDonald, Fensham, McAlpine, & Lawes, 2016), restoration thinning may dramatically accelerate the recovery of regrowth Brigalow forests to a state resembling mature forests (Dwyer et al., 2010a).
Many extensive regrowth ecosystems, including brigalow, span large climate gradients that extend into semi-arid zones. In such systems, the efficacy of restoration thinning is likely to depend strongly on the amount and variability of rainfall. For example, using a regional-scale survey of brigalow regrowth of different ages, Dwyer, Fensham, and Buckley (2010b) found that stand structural development and biomass accumulation were fastest in areas with high mean annual precipitation. However, it remains unclear exactly how precipitation modulates the ability of stems to respond to restoration thinning. Such knowledge will inform the spatial prioritisation of regrowth protection and management, especially in the context of ongoing climate change, which is anticipated to increase temperatures and reduce winter and spring rainfall in the region over the next 60 years (Ekström et al., 2015).
The present study capitalizes on an existing restoration thinning experiment established in 2007 in dense brigalow regrowth. The shortterm (0-2 years) growth and survival responses of the dominant A.
harpophylla were reported previously (Dwyer et al., 2010a); however, longer term data are required to reliably predict trajectories of biomass accumulation and structural recovery in response to potentially changing rainfall. In addition, co-occurring tree species such as C. cristata commonly persist in dense brigalow regrowth, albeit at much lower relative abundances than in original forests (Dwyer & Mason, 2018). Little is known about the responses of these once-codominant species to restoration thinning in dense brigalow regrowth. As such, we specifically address the following research questions: species, but C. cristata contributes the most above-ground biomass (Chandler, Buckley, & Dwyer, 2007). Prior to thinning the average stem density was ∼16,000 stems ha −1 compared to 1,500 stems ha −1 in the adjacent mature forest. Typical of brigalow regrowth, more than 99% of stems were A. harpophylla, and C.
cristata and was present only as scattered individuals that resprouted at low densities after the last clearing attempt (Dwyer & Mason, 2018).
Other tree species were absent or at extremely low abundance in the regrowth and were omitted as focal species for statistical analysis.

Experimental design
The experiment is a complete randomized block design, combined with pre-experimental sampling to permit Before-After-Control-Impact (BACI) inferences. In brief, sixteen 25 × 25 m plots nested within four blocks (i.e. four plots per block) were established in late 2006. These plots were large enough to capture an average of 1,000 stems before thinning. All plots were separated by at least 5 m and were >20 m from clearings. All woody stems of ≥1-cm diameter at 30 cm above ground level (AGL) were mapped and their circumference measured to the nearest centimetre. Due to the difficulty in determining connectivity between stems without excavation, stems were recorded as separate individuals if they were not visibly connected above the ground.

Statistical analyses
All statistical analyses were conducted in R Version 3.6.3 (R Core Team 2019) via RStudio Version 1.2.5033 (RStudio Team 2019). We fitted Bayesian models in Stan (Carpenter et al., 2017), via the brms package in R (Bürkner, 2017). In all cases, we adopted the default, weakly informative prior distributions for regression coefficients and variance parameters. As such, our models are analogous to generalized linear mixed effects models that are used widely in ecology.
We prefer Bayesian inference, however, because it allows for more flexible model fitting and generates posterior distributions for all model parameters. Significance of terms can be interpreted based on whether 95% credible intervals (CI) overlap zero. For all models, four Markov chains were used with a minimum of 2,000 iterations including a 1,000 iteration warm-up. Model convergence was assessed visually and via theR statistic (Gelman & Rubin, 1992). Explanatory variables were mean-centred (after transformation where necessary) in all models to facilitate convergence and aid interpretation of model coefficients.

Treatment effects on diameter growth from 2007 to 2017
Given the BACI design, we first tested for treatment differences in A. harpophylla diameter in both 2007 and 2017 using a hierarchical Bayesian model with lognormal errors. Although our primary interest was in testing for treatment differences that had emerged by 2017, we included 2007 diameters (i.e. before thinning occurred) to enable comparisons with initial stem sizes, and to ensure that pre-existing differences in diameters between treatments were negligible. Stem diameter was modelled as a function of treatment (four categories) and year (2007 vs. 2017). Random intercepts were assigned to each stem nested within each plot within each block to capture unexplained variation in stem, plot and block-level differences in diameter.
Because 2017 diameters exhibited visibly greater variance in thinned treatments compared to the control, we also allowed residual variation (σ) to vary by treatment for the 2017 measurements (all 2007 pre-thinning diameters were included in the control group for this purpose). Pair-wise differences between treatment groups in both 2007 and 2017 were assessed by calculating 95% CIs for each difference using the hypothesis function in brms.

Model of A. harpophylla growth
The growth of A. harpophylla between monitoring efforts was modelled as a function of the diameter of each stem at the start of each growth period (to capture size-dependent growth), precipitation during each period, neighbourhood stem density and two-and three-way interactions between these variables. Because there were multiple observations of growth for each stem through time, stem identity, plot and block were included as nested random intercepts. In addition, we fitted random slopes with respect to stem size (preceding diameter) for each stem. Finally, we allowed the residual standard deviation (σ) to vary by treatment given the visibly greater variation in growth increments observed in the thinned treatments compared to the control.
Growth was expressed as annual basal area increments (cm 2 year −1 ). In cases where a stem was recorded as dead in a given year, basal area increments were only calculated until the last measurement in which it was alive. Because stems were not measured every year, we divided the basal area increment between two measurements by the number of days in the increment period and multiplied by 365 to give the annual basal area growth rate. Basal area increments were negative in 2.6% of cases (55 of 2,101 increments) due to a combination of small measurement error and processes affecting stem diameters that are unrelated to growth (Chan et al., 2016). ELPD based on the widely applicable information criterion (WAIC) is a measure of the predictive accuracy of Bayesian models where models with higher ELPD are considered to be more accurate (Vehtari et al., 2017).

Model of A. harpophylla survival
The high proportion of surviving individuals at the end of the study (∼85%) prevented analysis of the effect of rainfall on the probability of survival through time. Thus, we modelled the probability that a stem was alive in 2017 as a function of its maximum diameter throughout the experiment and the neighbourhood density only using a Bayesian generalized linear mixed effects model with a Bernoulli distribution, including plot and block as random intercepts.

Effects of A. harpophylla density on growth rates of C. cristata
For this analysis, we quantified the relationship between A. harpophylla density and basal area growth of scattered C. cristata. We also quanti- harpophylla stems and their two-and three-way interactions. The main term of species identity allowed us to account for average differences in growth rate between species, whereas the three-way interaction allowed us to determine whether size-dependent growth rates of each species differentially responded to the density of neighbouring A. harpophylla stems. In addition to fitting random intercepts for each plot, nested within block, the effects of species identity, neighbourhood density and starting diameter were also permitted to vary by plot.

Treatment effects from 2007 to 2017
The diameter of A. harpophylla stems in 2017 was significantly higher in the 1,000, 2,000 and 4,000 stems ha −1 treatments than in the control plots (i.e. the 95% CI for the parameter estimates did not bound zero; Figure 1 & Table S1). Pairwise tests revealed that there were no significant differences between the 1,000 and ,2000 ha −1 treatments, but both had significantly larger diameters than the 4,000 stems ha −1 treatment in 2017. There were no significant differences in diameter between any of the treatments in 2007. In addition, in 2017, residual variation in diameter (σ) increased with the intensity of thinning.

Model of A. harpophylla growth
Of the 10 growth models tested, the neighbourhood density function with SD of 2 m and no weighting with respect to neighbour size had the lowest ELPD (Table S2). In this model, all main, two-way and threeway interaction coefficients were significant (Table 1). Neighbourhood density, rainfall and preceding diameter interacted in such a way that large stems in low-density neighbourhoods had the largest growth responses, and this effect was magnified in years with high rainfall ( When plotted on the scales of original units, the relationship between basal area accrual and preceding diameter was positive, but flattened off as diameter increased (Figures 2c and 2d), which was expected given the sqrt-transformed growth response and logtransformed diameter variable.
In terms of random effects, there was less unexplained variation among blocks and plots than among stems (both in terms of stemlevel intercepts and slopes). However, the most unexplained variation was within stems (i.e. departures from each stem's linear relationship with log(preceding diameter)). Predictably, this within-stem variation (σ) increased with the intensity of the thinning treatment (Table 1), reflecting greater unexplained variation in size-dependent growth in response to thinning than in unthinned, high-density controls.

Model of A. harpophylla survival
The probability of surviving to 2017 increased with the maximum diameter of stems (Table 2). In addition, there was a significant negative interaction between the maximum diameter of stems and neighbourhood density, such that smaller stems had a much lower probability of

Effects of A. harpophylla density on growth rates of C. cristata
For both C. cristata and A. harpophylla, larger stems accrued larger annual basal area increments. However, the significant three-way TA B L E 1 Summary of the Acacia harpophylla growth model including regression coefficients and random effect standard deviations

F I G U R E 2
Fitted relationships between annual basal area increments (cm 2 year −1 ) and the diameter of stems at the start of each growing period (preceding diameter [cm]) from the model of Acacia harpophylla growth, for both a low-rainfall year (a and c; 500 mm year −1 ) and a high-rainfall year (b and d; 700 mm year −1 ). Panels a and b show fitted relationships for the sqrt-transformed response variable and log-transformed diameter variable, as fitted in the model. Panels c and d show the relationships fitted on the original scales. In each plot, separate lines are fitted for each treatment using the average value of the sqrt(neighbourhood density) variable in each treatment. Points in panels a and c are raw data from the lowest third of annual rainfall values (< 501 mm year −1 ) and in panels b and d they are from the highest third of rainfall observation (> 611 mm year −1 ). Envelopes around each line are 95% credible intervals (most translucent) and 50% credible intervals (least translucent) calculated using only fixed effects (i.e. they are population-level credible intervals)

F I G U R E 3
Fitted relationships between the probability of survival and maximum diameter of each stem (cm) from the model of Acacia harpophylla survival. Separate lines are fitted for each treatment using the average value of the sqrt(neighbourhood density) variable in each treatment. Envelopes around each line are 95% credible intervals (most translucent) and 50% credible intervals (least translucent) calculated using only fixed effects (i.e. they are population-level credible intervals). To provide an indication of how the binary survival data correspond to each fitted relationship, we calculated points and standard error bars from nine bins of ordered binary data interaction between species identity, starting diameter and neighbourhood density (Table S3) revealed large differences between species in their growth-dependent responses to A. harpophylla density. For A.
harpophylla, the strong negative effect of conspecific density on sizedependent growth mirrored the growth model described above. For C.
cristata on the other hand, the slope of the relationship between diameter and growth rate was relatively insensitive to A. harpophylla density ( Figure 4). Instead, the overall density effect (Table S3)

DISCUSSION
This study combined an existing BACI restoration thinning experiment with temporal variation in rainfall to quantify the interactive effects of neighbourhood density and moisture availability on the growth and survival of A. harpophylla. It also compared the strength of A. harpophylla density effects on itself to its effects on C. cristata. Our results provide clear evidence that thinning permits A. harpophylla to grow rapidly during periods of high rainfall, and that interspecific competition between A. harpophylla and C. cristata is relatively weak.

4.1
What is the effect of the thinning treatment on A. harpophylla diameters a decade after thinning?
Reporting on the first 2 years of growth data from the present thinning trial, Dwyer et al. (2010a) found significantly higher growth in thinned plots compared to control plots. Ten years after thinning, the diameter of monitored A. harpophylla stems continued to be significantly larger in more intensely thinned plots. This confirms that accelerated stemlevel growth induced by restoration thinning continues over longer periods of time. In addition, we found that diameter variation increases with thinning intensity, conforming with similar observations made in a thinning trial in mixed species Eucalypt forest (Kariuki, 2008). This most likely reflects the wide range of neighbourhood densities created by random thinning, with some neighbourhoods remaining relatively dense and others becoming sparse.

How does neighbourhood stem density interact with interannual variation in precipitation to influence the growth and survival of retained A. harpophylla stems?
The kernel-smoothed stem density variables proved effective at explaining variation in A. harpophylla growth, providing further strong evidence that intraspecific competition is indeed limiting growth in dense brigalow stands (Dwyer et al., 2010a). The SD values used to generate the various neighbourhood density variables also provided insights about the scale of competition. Of the unweighted density variables examined, an SD of 1 m performed the worst by far, indicating that strong density effects extend beyond 1 m. An SD of 2 m performed best, suggesting that most density effects (i.e. 68.3%) occur within 2 m, although the model including an SD of 2.5 m performed similarly well (Table S2).
Thinning to 1,000 stems ha −1 (16% of pre-thinning basal area) produced very similar growth responses to the less severe 2,000 stems ha −1 treatment (26% of pre-thinning basal area). It could be that more shrubs have recruited into 1,000 stems ha −1 plots resulting in additional competition; however, a recent study of shrub recruitment in these plots found no difference in the number and size of recruits between the 1,000 and 2,000 stems ha −1 treatments (Dwyer & Mason, 2018). More likely, the similar growth responses indicate that the effects of intraspecific competition are minimal beyond a particular threshold of available space (Mainwaring & Maguire, 2004). As stems continue to grow, this threshold may increase and allow stems in the 1,000 stems ha −1 treatment to escape competition for longer than in the 2,000 stems ha −1 treatment, but this remains to be seen.
Predictably, the main effect of mean daily rainfall was positive and significant. However, consistent with other studies in forest ecosystems, the effect of rainfall on growth was strongly dependent on neighbourhood density, whereby stems experiencing less crowding were able to capitalize more strongly on periods of high rainfall (Dorman, Perevolotsky, Sarris, & Svoray, 2015, Ford et al., 2016, Sohn et al., 2013. Thus, thinning allows a typical stem to accrue more basal area F I G U R E 4 Relationships between annual basal area increments (cm 2 year −1 ) and the diameter of stems in 2007 (cm) from the model of Acacia harpophylla and Casuarina cristata growth. Panel a shows fitted relationships for Acacia harpophylla growth with separate lines fitted for high-density A. harpophylla neighbourhoods (75th percentile of density values) and low-density neighbourhoods (25th percentile of density values). Panel b shows the same for C. cristata growth. Relationships are fitted on the original scales of the growth response variable and the diameter explanatory variable. Points for high-density and low-density A. harpophylla neighbours are raw data from the highest and lowest thirds of sqrt(neighbourhood density) values, respectively. Envelopes around each line are 95% credible intervals (most translucent) and 50% credible intervals (least translucent) calculated using only fixed effects (i.e. they are population-level credible intervals) in a dry year than it would in a wet year if left unthinned. In a wet year, thinning allows stems to grow almost four times faster (compared to control plots in wet years).
Like growth, survival was positively size dependent, but thinning reduced the survival probability of smaller stems. This may be an effect of 'thinning shock' , a physiological response to a severe reduction in neighbourhood density possibly induced by heightened exposure of remaining stems to drought stress which has been shown to disproportionately affect smaller trees in a north American pine forest (Simonin, Kolb, Montes-Helu, & Koch, 2006). Alternatively, the removal of more productive stems from a physiologically integrated network may have reduced survival of smaller, more dependent stems (Liu, Liu, & Dong, 2016). Interestingly, although most stem mortality occurred in the first year after thinning, the effect of neighbourhood density only emerged after analysing the survival data over 10 years. Thus, in comparison to stem growth, which was shown to respond to thinning in just two years, thinning-induced stem mortality may be a process which manifests over longer periods of time.

4.3
Does thinning of A. harpophylla increase growth rates of sparsely regenerating C. cristata stems?
Consistent with the majority of field experiments in forests and other plant communities, high densities of A. harpophylla reduced the growth of conspecifics much more than heterospecifics (Adler et al., 2018).
Unfortunately, we were unable to quantify the strength of intraspecific competition for C. cristata or examine how rainfall modulates competition between A. harpophylla and C. cristata. However, the relatively low impact of A. harpophylla crowding on C. cristata growth suggests that, regardless of the limiting resource in this system, these species occupy complementary resource niches.

Implications for management
Our results have clear implications for the prioritisation and management of brigalow regrowth, especially for the purpose of sequestering carbon in the face of ongoing climate change. Under the medium emissions scenario (RCP 4.5), our study location is predicted to experience a 5% decline in annual rainfall and a 2 • C increase in temperature (CSIRO and the Australian Bureau of Meteorology 2019). Our growth model indicates that the effects on growth would be small in unthinned plots where competition inhibits growth even in wet years, but could be substantial in thinned plots (Figure 2). If these rainfall and temperature changes combined to exert moisture stress equivalent to a 50 mm year −1 reduction in annual rainfall, over 50 years a 7-cm-diameter stem in a thinned plot would grow to 21.5 cm instead of 26 cm diameter, resulting in 93 kg less above-ground biomass per stem (using Scanlan's [1991] allometric equation for A. harpophylla). This suggests substantial reductions in stand-level carbon sequestration under ongoing climate change.
Our findings of temporal responses to variations in rainfall likely apply to spatial rainfall gradients. With a mean annual rainfall of 620 mm, our study site lies in the middle of the regional rainfall gradient once dominated by brigalow forests and woodlands (Isbell 1962). Our results clearly indicate that thinning in higher rainfall areas will result in the greatest acceleration of growth. However, considerable ecotypic variation has been observed in A. harpophylla traits along rainfall gradients (Coaldrake, 1971, Johnson, 1964, which may also contribute to regional variation in growth rates and is worthy of further investigation. Our finding that C. cristata experiences relatively weak competition from A. harpophylla suggests that the low relative abundance of this species compared to nearby mature forests may be a consequence of dispersal limitation rather than competitive exclusion. Thus, active reintroduction of this once-codominant species may offer benefits for carbon sequestration through overyielding (Fichtner et al., 2018), if these species do indeed occupy complimentary niches. This could be achieved cost-effectively via direct seeding of C. cristata during highrainfall periods that are known to promote recruitment in natural populations (Chesterfield & Parsons, 1985), but further experimentation is required to determine how best to accelerate recruitment of formerly co-dominant tree species.

CONCLUSION
Our results show that randomly removing ∼75% of basal area from dense brigalow regrowth substantially accelerates the growth of retained stems, especially larger stems in wet years. This indicates that restoration thinning may be most effective at accelerating forest recovery in high-rainfall regions. But regardless of the location and thinning prescription, a safe and cost-effective alternative to manual ringbarking is required if thinning is to be implemented at meaningful scales.

AUTHORS' CONTRIBUTIONS
JMD conceived the idea and designed methodology. IRT and JMD collected the data. IRT and JMD analyzed the data. IRT and JMD led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.