Changes in liana density over 30 years in a Bornean rain forest supports the escape hypothesis

In two permanent plots of lowland dipterocarp forest at Danum a liana census in 1988 was repeated in 2018 using the same method. Trees ≥ 30 cm gbh (girth at breast height) were recorded for number of lianas ≥ 2 cm gbh on their stems. The forest was evidently in a late stage of recovery from a large natural disturbance. Over the interval, mean number of lianas per tree decreased by 22 and 34 % in plots 1 and 2. By 2018 there were relatively more trees with few lianas and relatively fewer trees with many lianas than in 1988. The redistribution was strongest for overstorey trees of the Dipterocarpaceae (more with no lianas by 2018) and understorey trees of the Euphorbiaceae (many losing high loads in plot 2). Proportion of trees with lianas rose overall by 3.5%. Number of lianas per tree showed a quadratic relationship with tree size (ln[gbh]): maximal for large trees, fewer for smaller and very large trees. Tree survival and stem growth rate were significantly negatively related to number of lianas, after accounting for spatial autocorrelation. Monte Carlo random samples of half of trees in 1988 were compared with the other half of trees in 2018. Relative frequency distributions differed significantly over time, but dipterocarps and euphorbs varied noticeably in their liana dynamics between plots. Regressions achieved best significant fits when number of lianas was a function of date, ln(gbh) and ln(gbh)2, but differently in the plots reflecting complicated host-liana dynamics. Analysis of most abundant trees species, individually, highlighted a group of emergent dipterocarps with low liana counts decreasing with time. Building on an earlier hypothesis, these trees lose their lianas with branch shedding, as they move into, and emerge from, the main canopy. They escape from the parasite. The process may in part explain the uneven nature of the forest canopy at Danum. Change in liana density was contingent on forest history and site succession, and plot-level structure and dynamics. Liana promotion in intermittent dry periods was seemingly being offset by closing of the forest and continued dominance by the Dipterocarpaceae.


INTRODUCTION
. Three climate-related processes or factors have been proposed to explain liana increases in neotropical, primary forests (i) Increased rates of natural perturbation and disturbance due to droughts and storms, which open up canopies and create more gaps; (ii) Increased evapotranspiration due to longer, more intense, dry seasons, 5 which favour lianas because of their deep rooting, stem anatomical adaptations, and highly efficient use of water, and (iii) increased global atmospheric carbon dioxide levels, which promotes the growth of lianas faster than it does of trees, especially under more lighted conditions; and (Schnitzer 2005, Swaine and Grace 2007, Cai et al. 2009, Schnitzer and Bongers 2011. 10 Lianas are nevertheless an important feature of aseasonal and less-recently disturbed rain forests in the paleotropics, particularly in SE Asia where primary forests are composed in recent decades largely of late-successional stages (Ashton 1964(Ashton , 2014Ashton and Hall 1992).
Liana abundance results here from the outcome of recent decades of structural stability combined with, in places, the legacy of major destabilizing natural events in the past (Newbery 15 et al. 1992(Newbery 15 et al. , 1999. When such a forest is in the late stages of recovery after a far-back major disturbance, a decline in liana abundance is expected (Campbell and Newbery 1993). The reasoning behind this postulate is that, as tree basal area abundance tends towards a site maximum, the main canopy closes, gaps become fewer and light conditions internal to the stand will be less favorable for liana recruitment and growth (Whitmore 1984). 20 With dry season seasonality and associated strong drought not being the issue, two factors remain though that could counteract an end-succession decline in SE Asian forests, namely a continuing rise in atmospheric carbon dioxide (Malhi and Phillips 2005) and perturbations caused by the intermittent effects of the El Nino Southern Oscillation, ENSO, (Walsh 1996b, Wright 2005b which opens the canopy, and leads the understorey temporarily 25 to have more lighted and drier conditions (Walsh andNewbery 2009, Newbery et al. 2011).
Lianas are also extensively distributed in logged forests of SE Asia (e.g. Pinard and Putz 1994), and have a major effect on forest regeneration; but that is a separate process from the dynamics of the mature primary forests, except for when these forests were last heavily disturbed by similar-scale natural events e.g. by extensive drought and/or fire (Beaman et al. 1985). Wright et al. (2015) reported on changes in liana abundance at Pasoh, Peninsular Malaysia, between 2002 and 2014, in terms of canopy cover infestation on trees ≥ 30 cm dbh 5 (95.4 cm gbh). They confirmed the lower susceptibility of dipterocarps to lianas compared with other tree families shown by Campbell and Newbery (1993). Tree mortality for all trees was affected when liana cover exceeded 75%. The change in proportions of all survivor trees infested was very small, declining from 52.3 to 47.9%. The frequency of trees with ≥ 50% cover of lianas barely altered. The samples were not independent though, because of the 10 surviving trees in common.
An analysis of the change in liana infestation over a 30-year period in a latesuccessional primary lowland dipterocarp forest in N.E. Borneo is presented here to test expectations that liana infestation decreased over time and lianas were still having negative effects on tree growth and survival. The hypothesis of Campbell and Newbery (1993), that liana 15 loss is largely controlled by branch fall for overstorey trees and whole tree death for understorey ones, is extended to show that this mechanism allows especially the dipterocarps to escape their parasite load when reaching the upper canopy and becoming emergent.
with an average monthly temperature of 26.9°C and an average monthly rainfall of 2832 mm (1986( -2007( , Newbery et al. 2011. More detailed information is given in Newbery et al. (1992Newbery et al. ( , 1996 and Walsh and Newbery (1999). The two plots are designated 'main' plots, MP1 and MP2, at Danum to distinguish them from 10 smaller neighbouring satellite ones. 5 The permanent plots were set up between 1985 and 1986 and have since been under ongoing long-term observation (Newbery et al. 1992). All trees with a girth at breast height ≥ 10 cm (gbh, at 1.3 m or above buttresses) were tagged, mapped, identified and gbh was measured. Since then, four more censuses have been made, the latest in 2015 (D. M. Newbery et al., unpubl. data). In each of those censuses, survivors, recruits and dead trees 10 were recorded and gbh measured for all live trees. The two plots were used as replicates for the present study, since they were shown to be similar in species richness and forest structure (Newbery et al. 1992). This means that analyses were performed, unless specified differently, for each plot separately. From September 1987to February 1988and June to September 1989(mid-date September 1988, Campbell and Newbery (1993) conducted a liana census in the two plots.

Fieldwork
Liana stems on all trees ≥ 30 cm gbh were counted (Zahnd 2018). The present study in June and July 2018, 30 years later, closely followed the procedure of the first one to achieve the 20 highest feasible level of compatibility. Liana stems ≥ 2 cm gbh on all trees ≥ 30 cm gbh were again counted. Lianas were defined here as woody climbing plants which start growing from the ground. Lianas were regarded lianas as individuals when they rooted close to the tree observed, or when crossing over into the target tree from another tree. Only lianas which had foliage in the canopy of the tree were counted: lianas which did not reach the crown of the tree, or only made a loop and descended, or that were clearly dead, were not counted. As in the 1988 census, lianas were not identified taxonomically because taking botanical samples for identification would have influenced the growth and mortality of the trees themselves undergoing long-term observation (Campbell and Newbery 1993). 5 For the present study tree girths measured in the first (1986) and last (2015) main plot tree censuses were applied to the liana counts from 1988 and 2018. These censuses gave the tree data closest in terms of plot census dates to those of the liana censuses. For both liana censuses, trees which died between 1986 and 1988 or between 2015 and 2018 respectively were disregarded, as were trees which only grew bigger than 30 cm gbh in the respective size structure, forest dynamics and species composition, even in moderately small ways, combining them into one sample would have confounded several important and interesting differences that bear on the relationships between liana density and dynamics and these factors. Moreover, the plots are replicates at the 4-ha scale and separated enough spatially to 5 provide some indication of forest variability on the local landscape. Apart from revisions to nomenclature, the 1988 liana data used here are the same as those reported in Campbell and Newbery (1993).
Trees at each census were divided into four size classes: 1, 30 -< 60; 2, 60 -< 120; 3, 120 -< 240; and 4, ≥ 240 cm gbh. These are referred to as medium-sized, large, very large 10 and largest trees respectively: the class 'small' would be those 10 -< 30 cm gbh, not considered for the liana censuses. Over-understory index (OUI, scale 0 -100), for the 100 most abundant species in the plots, was adopted from Newbery et al. (2011). It is based on first axis of a principal components analyses of ratios of numbers of trees ≥ 30 cm gbh/≥ 10 cm gbh and ratios of basal areas of the same size classes, the rescaled scores averaged over 15 four plot tree censuses 1986 to 2007. Three storeys were designated: overstorey (OUI > 55), intermediate  and understorey (OUI < 20). The 30-cm-gbh threshold is the same as the lower one for liana recording.
Analyses of the liana census data were mainly done in the computing language R (R_Core_Team. 2020;and Fox and Weisberg 2011). Special commands are listed in 20 Appendix 1, where part 1 has the basic code for randomization testing of differences over time, and part 2 lays out the essential regression models, particularly those that allowed for spatial autocorrelation. 25

Error distributions
The negative binomial distribution fitted the frequencies of number of lianas per tree far better than the Poisson, with χ 2 -values of 19.34 and 14.41 in 1988 for MP1 and MP2 respectively (df = 7), and correspondingly 13.89 and 14.50 in 2018 (df = 5). All four were 5 marginally significant at P ≤ 0.05. For comparison, fitting the Poisson χ 2 -values were 3344 and 1185 in 1988 and 2018 (plots together). The k-values (expressing degree of aggregation when < 1) were, for the four cases in order: 0.625, 0.491, 1.265 and 1.034; indicating that lianas were much more aggregated in 1988 than in 2018. The geometric Poisson (Polya-Aeppli) distribution gave very similar fits to the negative binomial, and since the latter is 10 fairly widely incorporated as a family error distribution in most regression software, it was taken the error model. The remaining differences in fit were largely due to observed frequencies being less than expected ones for counts of 1 and 2 and the converse for counts of 3 and 4, for which no single-mode function can readily cater.

Numbers of lianas
Negative binomial GLM regressions were made using the 'glm.nb' function in MASS library of R). Comparing fits (by anova) of nr_lianas ~ ln(gbh) and nr_lianas ~ ln(gbh) + elev for 1988 led to LR's (log-likelihood ratios) of 11.86 (P ≤ 0.001) and 9.32 (P ≤ 0.01) in MP1 and MP2 respectively, and correspondingly for 2018 LR's of 13.34 and 40.43 (both P ≤ 20 0.001), indicating that including elevation enhanced the model significantly. Here and later in model formulations, the term 'nr_lianas' is number of lianas per tree. Adding a ln(gbh)•elev interaction to the two-term model led to no improvement, LR's of 0.14 and 0.49 in 1988 and 0.02 and < 0.01 in 2018 (all P > 0.10). By contrast, including instead a [ln(gbh)] 2 term led in 1988 to LR's of 3.10 (P ≤ 0.1) and 31.87 (P ≤ 0.001) in MP1 and MP2, and 2018 LR's of striations for discrete data: residuals themselves however were positively skewed.
To test and cater for spatial autocorrelation h-likelihood GLM (HGLM) was run using the 'hglm' function in same-named hglm library of R, by including a SAR correlation-term based on plot inverse-distance matrices. Hglm does not cater for a negative binomial error so this had to be substituted by a quasi-poisson one: and models without the spatial 10 autocorrelation adjustments compared with the negative binomial one. Running the selected model with the interaction also resulted in later being insignificant for all census x plot combinations. The method is described in Alam et al. (2015), Lee et al. (2017a, b) and Rönnegard et al. (2010).

Tree survival
Tree survival between 1988 and 2018 was modelled with a logistic GLM regression (Hilbe 2009(Hilbe , 2011, as a standard 'glm' (binomial error) starting with survival ~ sqrt(nr_lianas88). The transformation normalized the original nr_lianas variable to a large extent but not entirely. Adding ln(gbh88) and elevation resulted in varying degrees of 20 improvement to the fits depending on the plot. With ln(gbh) the LR's were 4.50 (P ≤ 0.05) and 2.67 (P > 0.10) in MP1 and MP2 respectively, whilst adding elevation led correspondingly to LR's of 0.02 (P > 0.05) and 7.47 (P ≤ 0.01), although including both terms they were 4.58 (P ≤ 0.05) and 10.37 (P ≤ 0.01). As a compromise, in order to compare plots on a same-model basis, all three terms were retained to have a final model: survival ~ sqrt(nr_lianas88) + ln(gbh88) + elev.
To account for spatial autocorrelation, an autologistic regression was run with the command 'logistic.regression' in the library spatialEco in R, using an inverse squared-5 distance matrix as a measure of correlation. Further, as an alternative to handling spatial autocorrelation affecting the model, repeated stratified sampling was achieved by taking one tree at random per 20-m x 20-m sub-plot (100 trees per plot), N' = 500 times, and averaging the coefficients and statistics. This meant that trees on average 20-m apart were assumed to have independent liana counts. The repeated runs will have reused trees many times and not 10 be independent, especially when there were very few per sub-plot from which to sample. This additional approach was intended to be confirmatory of the regressions, free of distancedecay function. Methods are outlined and discussed in Besag (1972), Cliff & Ord (1981), Augustin et al. (1996), Dormann (2007 and Zuur et al. (2009). Analysis at the family level, for the Dipterocarpaceae and Euphorbiaceae, failed to 15 converge to solution with the autologistic function in three of four cases, and it was replaced by a GLMM using the coordinates of 20-m x 20-m subplots as the cluster term in command 'glmmML' of the library glmmML in R (Rhodes et al. 2009, Brostrom andHolmberg 2011).

Tree growth 20
With standard linear, or equivalently generalized least squares (GLS), regression (gaussian error) as the function 'lm' in R, relative growth rates (rgr, mm/m/yr) of trees, that is survivors, between 1988 and 2018, were fitted first as ~ sqrt(nr_lianas88). Adding the term ln(gbh88) considerably improved model fits, with change in F-ratio (equivalent to LR) of 39.81 and 30.51 (both P ≤ 0.001) in MP1 and MP2 respectively, but not consistently so much in the case of elevation, the corresponding LR's being 3.60 (P ≤ 0.1) and 0.62 (P > 0.10): both terms included led to combined improvements over sqrt(nr_lianas88) alone, LR's of 21.08 and 15.53 (both P ≤ 0.001). Adding the squared size-term, [ln(gbh88)] 2 , to ln(gbh88) and sqrt(nr_liana88) was of little gain, with F-ratios of 0.10 and 2.10, and adding to these 5 elevation, 1.20 and 1.28 (all four at P > 0.10). Attempting to cater for some remaining nonnormality in rgr, by using gamma or inverse-gaussian errors, resulted in scarce improvement in either fits or residual plots. Accordingly, the best parsimonious model was: rgr ~ sqrt(nr_lianas88) + ln(gbh88).
With a rgr regression using the gaussian error, there were more ways of dealing with 10 spatial autocorrelation than those (currently as software) available for regressions with binomial and negative binomial errors. Firstly, spatial dependence was modelled using a variogram of exponential form in GLS (a correlation structure), within the function 'gls' in library lnme in R. Comparing the accepted model form, without and with the spatial correlation term, led to modest improvement in the fits, with LR's of 8.84 (P ≤ 0.05) and 5.54 15 (P ≤ 0.10) for MP1 and MP2 respectively. Using a spherical form in the variogram made hardly any difference. Randomization runs were conducted in the same way as they were for the survival data. A different, additional, approach was to use a correlation matrix again based on inverse distances but guided by prior analysis of spatial autocorrelation via a correlogram. Two indices, Moran's I and Geary's C, were calculated at 5-m interval distances 20 between 0 and 40 m. Moran's I for MP1 was significant at 5-10 m (P ≤ 0.001) and at 30-35 m (P ≤ 0.1), and for MP2 at 0-5 m (P ≤ 0.1), 10-15 and 15-20 m (both P ≤ 0.05). Geary's C was significant in MP1 at 0-5 m (P ≤ 0.01) and 5-10 m too (P ≤ 0.05), and in MP2 again at 0-5 m (P ≤ 0.001) and 15-20 m (P ≤ 0.10). Thus, up to 20 m, both indices show significant spatial correlation, permitting use of the 'spautolm' function in library spatialreg in R, with SAR mode and range 0-20 m. The methods used here follow Pinheiro and Bates (2000), Haining (2010), Plant (2012), Bivand et al. (2013) and Borcard et al. (2018).

Change in liana abundance with time 5
Changes in the proportion of trees carrying lianas, the frequencies of trees in liana density classes, and the mean number of lianas per tree, between 1988 and 2018 were estimated. Some of the trees were different ones at the two dates (those which died after 1988 were lost and those which recruited after 1988 to survive to 2018 were gained), but close to half were the same individuals (the survivors). (The 1721 survivors formed 51.2% of the trees 10 in 1988 and 48.7% of those in 2018.) Therefore, the tree populations at the two dates were not fully statistically independent. Subsamples of different trees were accordingly selected for 1988 and for 2018, by taking at random half of those in 1988, flagging them, and the unflagged trees contributed to the sample for 2018. Some unflagged trees will not have survived 1988 to 2018, and they were replaced (approximately in terms of numbers) by 15 recruits appearing in 2018. Thus no one tree was present in both dates' subsamples. Statistical tests were repeated N' = 500 times as a Monte Carlo procedure (Manly 1997). This randomization procedure not only eliminated temporal dependence but it also reduced the effects of spatial autocorrelation.
Differences in the proportions of trees having lianas, that is presence or absence of 20 lianas, at the two dates were tested using Pearson's chi-squared statistic (2 x 2 test; df = 1).
For frequencies of trees in liana density classes (0, 1, 2…r lianas) the chi-squared statistic was again applied (2 x r test, df = r − 1), with the provision, achieved by pooling higher tail classes, that no class had < 5 trees. Firstly, all trees were used in single tests, without any subsampling. Randomization tests were then run N' times, the mean statistic found and its significance, and the null hypothesis was rejected overall when ≥ 80 % of the randomization runs were individually significant at the P ≤ 0.05. The mean chi-squared values were taken as representing 'typical' differences in proportions free of any temporal autocorrelation effects.
Additionally, McNemar's chi-squared test of symmetry (Agresti 2007), was employed to ask 5 whether proportions of trees with or without lianas changed with date, for just the 1988survivors. The symmetry statistic could only be calculated for counts in the density classes when taking all trees: the two main families had insufficient numbers.
Change in proportions of trees with lianas between the two dates could be analysed too with the negative binomial GLM, to ask whether mean numbers per tree differed 10 significantly, and whether the relationship between proportion and size of tree changed over time. Regressions were run with subsamples at the two dates from the same N' = 500 randomizations; they compared nested models 1 to 4: nr_lianas ~ year, ~ year + ln(gbh), ~ year * ln(gbh), and ~ year * ln(gbh) + [ln(gbh)] 2 + year: [ln(gbh)] 2 . The last two models tested whether differences in proportions were further accountable for by linear and non-15 linear interactions between year and ln(gbh). Interpreting the interactions in GLMs between a categorical factor (e.g. year) and continuous variables (ln [gbh] and [ln(gbh)] 2 ), and then using the term to infer estimates of the expected response variable (nr_lianas) requires special consideration. Guidance to the calculations are found in Hilbe (2009), Cameron and Trivedi (2013) and Long and Freese (2014).

Liana abundance and dynamic status
Starting with the trees censused in 1988, proportions of trees with lianas could be compared between those dying and those surviving by 2018 with the chi-squared test of association, and the numbers of lianas per tree similarly compared with negative binomial 25 GLM. In complementary way, those trees at 2018 that were survivors from 1988 and those that had recruited (into the ≥ 30 cm gbh population) by then could also be compared. for period one, between the first two tree censuses (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996), have been reported in

Liana frequency distributions
The nine to ten most abundant families (Table 1 and Appendix 2: Table S1) combined 20 made up 75% and 71% of all trees in MP1 and MP2 respectively. The most abundant families were the Dipterocarpaceae and the Euphorbiaceae, followed by the Meliaceae, Lauraceae, Phyllanthaceae and Annonaceae. Relative abundances of the first 13 families per plot plus all other pooled into a 14 th , did not change significantly between dates (MP1: χ 2 = 9.57, df = 13, P = 0.73; plot 2: χ 2 = 18.05, df = 13, P = 0.16), nor did the overall 10-cm gbh classed distributions of all trees (pooling trees ≥ 150 cm) between years (MP1: χ 2 = 14.19, df = 12, P = 0.286; MP2: χ 2 = 9.35, df = 12, P = 0.67; see also Appendix 2: Fig. S1). This largely rules out that changes in overall liana density were due to major shifts in the tree community over 5 the 30 years. Nevertheless, MP1 and MP2 did differ in one important respect: MP1 had more very large trees than MP2 (Appendix 2: Fig. S1), although they decreased in number more in MP2 than MP1 between 1988 and 2018. These very large trees were mostly dipterocarps.
In total, 2865 and 2442 liana stems were counted in plots 1 and 2, respectively. This is a 21.8% and 27.0% decrease since the first census (3688 and 3347 lianas in plots 1 and 2, 10 respectively). The mean number of lianas per tree, excluding trees without any lianas, decreased in MP1 and MP2 from respectively 2.16 and 2.02 in the first census to 1.69 and 1.33 in the second one, decreases of 21.8 and 34.2% (Table 1). For all trees, the proportion of liana-infested trees significantly increased by 4.7% in MP1 and by 2.4% in MP2. None of the major families individually showed significant changes in the proportions, although in MP1 15 trees in the Dipterocarpaceae decreased whilst those in the Euphorbiaceae correspondingly increased: differences for these two families in MP2 were negligible (Table 1).
Frequency distributions of counts of lianas per tree declined smoothly in plots MP1 and MP2, with relatively more trees with no or few (1 to 3) lianas, and relatively fewer trees with many (≥ 7) lianas in 2018 than in 1988 (Fig. 1a, b). The distributions appeared 20 approximately exponential in form. For the Dipterocarpaceae inflation of zero counts was obvious, with the difference between years matching the overall one in MP2 but not so well in MP1 (Fig. 1c, d). The Euphorbiaceae followed a scaled down representation of the overall plot distributions (also for the zero-class), but MP1 not (Fig. 1e, f). The two plots appeared on first examination to be showing different dynamics in liana densities over time, which same as those used in Campbell and Newbery (1993: Figs 2 and 3), albeit presented here slightly differently. In that earlier paper the size distributions of trees in the Dipterocarpaceae and Euphorbiaceae were compared (ibid.: Fig. 1).
Considering the frequency distributions of numbers of lianas per tree in the four size 5 classes (Fig. 2), the gradual decline for 0 to ≥ 8 lianas per tree, in both plots and at both dates, in the medium-sized trees (class 1) becomes increasingly more discontinuous as proportionally more trees up to the very large trees (class 3) have no lianas, and among the largest (class 4), zero counts become dominant. This implies indirectly that as the trees become larger they are shedding their lianas.

Number of lianas and tree size
The relationship between number of lianas per tree and tree size (gbh) could be described by a quadratic form, for MP1 in 1988 rather weakly, but more strongly in 2018 (Table 2, Fig. 3). For MP2, the fits were similarly strongly quadratic at both census dates.

15
More lianas were recorded on medium than on small and large trees: the distinction between censuses in MP1 appeared mainly due to large trees in 1988 having similar loads to mediumsized ones. Consequently, change in numbers over time was more marked for MP1 than for MP2. In 1988, coefficients were statistically weaker for MP1 than for MP2; in 2018 they were highly significant for both plots. Adjustment for spatial autocorrelation was small 20 although its effect was more noticeable for MP1 in 1988, especially on the gbh-terms.
Regression using HGLM with the quasi-Poisson error gave similar results to the negative binomial GLM (Appendix 3: Tables S1a, b), suggesting that the variance-mean relationship in numbers of lianas was close to linear (see van Hoef and Boveng 2007).
Pseudo-R 2 -values (as %), based on differences in deviances for fitted and null models (in the latter the mean was the sole term), and their random and fixed components (in parenthesis), were in 1988, 24.2 (22.0, 2.2) and 16.8 (14.8, 2.0) for MP1 and MP2 respectively, and correspondingly in 2018, 20 (17.4 and 2.6) and 10.2 (4.6, 5.6). The  Table S1c). Accounting for spatial autocorrelation (HGLM) affected the coefficients only very slightly but the significances more, being lower for the Euphorbiaceae (remaining at P ≤ 0.05). Elevation was otherwise also significant for two other combinations (P ≤ 0.05). Trees at the family levels were less dense than overall, so the influence of spatial autocorrelation was less important.

Tree survival
Tree survival was significantly negatively related to number of lianas in both plots, although weakly with gbh and inconsistently (positively) with elevation in MP2 (Table 3).
Including a spatial correlation term (in the autologistic regression) changed the slope of survival on number of lianas very little, even though the spatial covariate was significant for MP1 (Appendix 3: Table S2a). The coefficient for gbh was likewise minimally affected.
Stratified random sampling gave average slopes slightly smaller than those for the autologistic fit, especially for the number of lianas term in MP1 (Appendix 3: Table S2b).

5
Since number of lianas in the regression were those recorded in 1988, it is reasonable to conclude that increasing liana loads reduced tree survival in the 30 years. Despite the clear trend, the pseudo-R 2 values were very low, 0.8 and 1.5% in MP1 and MP2 respectively. At the family level, survival was significantly negatively dependent on number of lianas for euphorbs in both plots (P ≤ 0.05) and dipterocarps in MP2 (P ≤ 0.01), although for MP1 the 10 slope was in the same direction (Appendix 3: Tables S2c, d). Interestingly, the slopes of all four relationships are higher than those for all species in the plots (Table 3). What the liana loads were later after 1988 on trees that died by 2018 is unknown.

15
The negative effects of number of lianas on tree relative growth rate (rgr88-18) were highly significant (P ≤ 0.001). Elevation was non-significant, but gbh more significant (also negative) than number of lianas (Table 4). Coefficients were very similar for MP1 and MP2.
Accounting for spatial autocorrelation led to scarcely any change in coefficient values, for either the two regression approaches or for stratified random sampling, indicating that local 20 spatial dependencies were very weak (Appendix 3: Tables S3a-c). Nevertheless, testing for autocorrelation avoided type I errors. The pseudo-R 2 -values for growth models were much higher than for survival, 22.5 and 32.0% for MP1 and MP2.
Regressing the number of lianas in 1988, with a fuller model that included ln(gbh) and [ln(gbh)] 2 of 1988, and elevation, the coefficient for rgr decreased for MP1 and increased 5 moderately for MP2 (Appendix 3: Table S4). A similar regression with gbh of 2018 instead, led to the rgr coefficients increasing much more for both plots, so that by partialling out the effects of gbh in the negative dependence on liana load on rgr increased. Lastly, the change in number of lianas between 1988 and 2018, treated as a Poisson variable (n + 22 -to make all changes positive) was insignificantly related to rgr in MP1 (−0.006 ± 0.028, z = 0.234, P 10 = 0.82) and in MP2 (0.024 ± 0.027, z = 0.879, P = 0.38). Adding ln(gbh) and elevation terms did not improve the fitting overall or significance of rgr.
Linear regression of rgr on number of lianas and ln(gbh) at the family level showed the expected dependence on size for the dipterocarps but not the euphorbs (Appendix 3: Table S3d, e) but largely no or weak significance of dependence on number of lianas. The use 15 of the variogram-based spatial autocorrelation adjustment scarcely affected the results: with coefficients increasing slightly, and for euphorbs in MP2 a small rise in significance (still only P ≤ 0.05).

Changes in liana abundance -contingency analysis 20
Using all trees at both census dates, there was a small significant (P < 0.01) increase of almost 5% in the proportion of trees with lianas between 1988 and 2018 for MP1 but not MP2 where it was less at ca. 2.5% (Table 5). At the family level there were no clear significant differences, largely because sample sizes were much smaller than for all plot trees.
The randomization testing supported the result for plot 1 but only at P < 0.05, and 80 percentile values were all insignificant (Table 5). Here a qualification is needed: χ 2 values were positively skewed, more so in MP2 than MP1, so that the medians were lower than the means. This made the significance for all trees in MP1 possibly more marginal. Considering just the survivors, the change in MP1 and also in MP2 was larger with 8-9% increases over 5 time, and more strongly significant (P < 0.001) than when deaths after 1988 and recruits at 2018 were involved. Between families, there was pronounced difference in MP1: with a strong decrease in proportion of trees in the Dipterocarpaceae (not quite significant however) counterbalanced by a much larger increase for the Euphorbiaceae (~22%). These differences roughly mirror the smaller relative changes at the 'all trees' level. However, for MP2 10 dipterocarps did not alter proportions and euphorbs increased just slightly (Table 6).
Therefore, over all trees the proportion with lianas increased slightly in the 30 years, more so for just survivors, but mostly for surviving euphorbs in MP1.
Considering counts of lianas per tree is more informative. Taking all trees (i.e. all families), and the single χ 2 -tests with no randomized subsampling, very highly significant 15 differences (P < 0.001) in the relative frequency distributions between dates are highlighted for both plots (Table 6, Fig. 1). The large χ 2values in general were close to normally distributed, and so the medians were also almost the same as the means. Not only do the randomization tests support the simpler ones, but also so do the McNemar's symmetry tests for just survivors. Figure 1 shows clearly the change in class frequencies with a decline in 20 numbers of trees with very high loads (≥ 8 lianas per tree), and an increase in trees with 1-3 lianas per tree. Considering the families separately, the randomization tests did not hold up the confidence in the single overall χ 2 -test (the latter significant at P < 0.05 in both plots, Table 6) for the dipterocarps. Conversely, whilst the Euphorbiaceae showed no evidence of differing class proportions in MP1, they did in MP2, supported by the randomization tests (P The failure to detect significance in the dipterocarps under the randomization may have been because there were few trees in very-large tree size class, the one with the high liana loads especially in 1988 and hence sample variances increased when taking only half of the trees at each census. Fig. 1 15 Negative binomial GLM regressions using the randomized subsamples involved four models, and for all trees each of them is shown in Table 7a. Model 3, with three terms (aside from the intercept), is interesting because, for MP1, year, size and their interaction were all significant. Without the interaction term, year became more significant but size less so, though more noticeable was that between models 2 and 3 the sign of the coefficient for year 20 changed from negative to positive and accordingly trends with ln(gbh) were catered for more in the latter by the interaction term. For MP2 it was just year that had a consistent and strong (negative) effect. Model 3 was expanded to model 4 to have two further terms, [ln(gbh)] 2 and its interaction with year, to accommodate the quadratic element of the relationship nr_lianas and ln(gbh) and thus the different shapes of the curves between years particularly for MP1.

Changes in liana abundance -regression analysis
Even with more terms in model 4, than the other models, individual terms were all less significant for MP1 but for MP2 ln(gbh) and [ln(gbh)] 2 featured strongly (Table 7a). It is important to recall that model 4 was more a descriptive than hypothesis-testing model.
Model 4 was overall the best-fitting model for MP1 and MP2 but for slightly differing 5 reasons. Averaged over the 500 iterations, of all four models, AIC was lowest in model 4 for both plots, and model 4 reached highest pseudo-R 2 values in MP1 and MP2 (2.64 and 4.56 % respectively) (Appendix 4: Table S1). Highest LRs were attained when model 4 was compared to model 1, again in both plots (LR03 = 32.3 and 46.4); the changes from model 1 to 3 were nearly as significant (LR13 = 22.7 and 45.8). The six paired comparisons were the 10 models coded 12, 13, 14, 23, 24 and 34. Notable was that LRs for models 1 to 2 in both plots were relatively low (Appendix 4: Table S1), indicating that model fitting markedly improved when interactions between year and ln(gbh) and [ln(gbh)] 2 were estimated. Probabilities of the χ 2 change in deviance when comparing models were lowest for model 1 versus 4 in MP1, and even lower for models 1, 2 and 3 with model 4 in MP2.

15
Considering the frequencies with which models were selected across the iterations, lowest AIC-values were attained with model 4 in 497 of 500 iterations for MP1 (model 3 in three), and all 500 for MP2. Likewise, model 4 had the highest LR-values among all six comparisons in all 500 for both plots, and similarly pseudo-R 2 was highest for all 500 per plot. Such consistency is not surprising given the small SEs of the statistics (Appendix 4: 20 Table S1). Delta-AIC was found for models paired in the same way done for LR, and here 475 of 500 were largest for model 1 minus model 4 for MP1, but 373 of 500 for model 3 minus model 4 in MP2 (Appendix 4: Table S2). Similarly, the significance of the χ 2 measure of deviance change was highest for 455 or 500 iterations in MP1 for the 14 comparison and 438 or 500 in MP2 for the 34 one. LR was never highest for the model 23 comparison or these statistics the highest respectively lowest for model comparisons 12 and 13 either (Appendix 4: Table S2). These last results point to an important difference between the plots: whilst model 4 was overall best for both, the most improvement for MP1 was between model 1 and 4, but for MP2 it was between models 3 and 4.

5
In these glm.nb fits, year was coded as a factor. This meant that year18 was the difference from the reference year88, and likewise interactions between ln(gbh) and [ln(gbh)] 2 with year were differences in their corresponding slopes (Appendix 3: Table S4). For MP1 nr_lianas in 2018 became a decreasing proportion of the nr_lianas in 1988 as gbh increased, going from 85% at ln(gbh) = 3.5 (gbh = 33.1 cm) to 40% at ln(gbh) = 5.5 (gbh = 244.7 cm) 10 (Appendix2: Table S5a). But for MP2 the change was very different, from 62 rising to 73% for the same gbh range (Appendix 3: Table S5b). So whilst model 4 was the better one for MP1 (cf. models 1-3) and showed the strongest change with respect to gbh, its main terms became less, and the interaction terms conversely more, significant.
Model 4 was also run on the full data, i.e. all trees without random subsampling at the 15 dates, and the coefficients fitted matched the averages from the randomizations very well.
(Appendix 3: Table S6a). Further, the coefficients from model 4 could be used to reconstruct the individual equations for the dependence of nr_lianas on ln(gbh) and [ln(gbh)] 2 at each date (similar to those reported in Table 2 without the elevation term, again matching closely (Appendix 3: Table S6b). This demonstrated that in testing for change in nr_lianas over time 20 (the year factor) the gbh curves were being modelled by the interaction terms correctly.
Predicted (marginal) means of nr_lianas in 1988, and in 2018 -using a mean ln(gbh) and [ln(gbh)] 2 common to years -were very similar across the four models, with just small decreases of 1.5-1.8% for MP1, and 2.8-3.0% in MP2, between models 1 and 4 (Appendix 3: Table S7a). Predicted means at the family level were also very similar (Appendix 3: Table  7b). The means resulting from models 1 and 2 are, naturally, virtually the same as the values given in Table 1.
Liana abundance and dynamic status 5 Among the trees that died 1988-2018, a higher proportion had lianas present in 1988 compared with those that survived (Table 8a). This was especially significant in MP2.
Correspondingly, dead trees had on average more lianas per tree than survivors. Trees in both the Dipterocarpaceae and Euphorbiaceae followed the same trends, again more strongly in MP2 than MP1. Survivors of dipterocarps had less than half the numbers per tree than those 10 that died. Recruits on the other hand by 2018 had proportionally fewer trees with lianas than survivors, and less lianas per tree, especially in MP1 (Table 8b). Differences were less strong at the family level, except for dipterocarps where the proportion of trees with lianas was significantly less on survivors that on recruits in both plots, though numbers per tree were similar.

Liana densities per tree for individual species
The mean numbers of lianas per tree for the 51 and 52 more frequent tree species in MP1 and MP2 together at 1988 and 2018 were merged to 61 species in common (Appendix 5: Table S1). Mean species' number of lianas per tree was very weakly correlated with OUI 20 (r = −0.006, df = 49, P = 0.97) and ln(gbh) (r = −0.047, df = 49, P = 0.74) in 1988 though slightly more strongly in 2018 (r = −0.175, df = 50, P = 0.21; r = −0.278, df = 50, P = 0.046 correspondingly). OUI and ln(gbh) were highly correlated in 1988 (r = 0.859, df = 48, P < 0.001; with Nothaphoebe sp. dropped because its ln(gbh) was very heavily outlying) and in 2018 (r = 0.904, df = 50, P < 0.001). However, the poor correlations hide an interesting separation between six to seven canopy/emergent dipterocarps, with high OUI-values but very low numbers of lianas per tree, versus the rest that follow a general positive trend (Fig.   4a, b). The separation is weak in 1988 but become distinct in 2018. Similar patterns arise when plotting against ln(gbh) instead of OUI. This separated group of dipterocarps clearly 5 have very much lower numbers of lianas per tree than would be predicted from the general trend with increasing OUI of all the others species.

Changing liana abundance
The two liana censuses covered a recent 30-year period in a forest considered to be in a stage of late successional recovery from a major natural drought disturbance, probably in 15 the 1870's. Indicators were the very small proportion of plot area with gaps, very few pioneer-stage trees, few large lying stems, no evidence of major wind damage, and total tree basal area approaching an expected site maximum accompanied by declining tree density (Newbery et al. 1992). This contrasts strongly with the majority of liana studies to date which were conducted in heavily disturbed and secondary successional forests, or forests with a  Large-sized trees had mainly higher liana loads than medium-sized and very large ones, except for MP1 in 1988 when very large trees had similarly high loads to large trees. The statistical errors were relatively high here due to the fewer trees in the very large-tree size class. Accordingly, this latter class showed a much faster decline between the censuses in 20 MP1 than in MP2. Number of lianas per tree decreased significantly with topographic elevation overall: they were fewer on ridges than lower slopes. Despite the often high significance of coefficients from the model fitting, total variances accounted for were low, indicating that many other unknown sources of variation besides tree size and elevation were important.
In both plots trees with many lianas in 1988 survived significantly less well than those with none or few. Tree survival also decreased with elevation, more strongly in MP2 than MP1, but it was unrelated to tree size. Likewise, when liana loads increased this led to significantly reduced rgr of gbh, but rgr was unaffected by elevation. Final model fits 5 accounted for far more variance in growth (27% on average) than in survival (~ 1%). If lianas were affecting tree growth negatively, then survival would be expected to decrease as well, as reduced growth normally enhances tree mortality. Partialling out elevation still left survival and growth depending on tree size. Differences between families in the dependence of their trees' growth rates on liana load were not consistent between plots though. For dipterocarps 10 the relationship was negative in MP2 but not so in MP1, and for euphorbs it was slightly negative both in MP1 and in MP2. The decrease in rgr with increasing tree size for all trees was highly significantly, for the dipterocarps too, but not for the euphorbs, at the family level.
Monte Carlo randomization testing showed that the increase in proportions of trees with lianas was marginally significant in MP1 but not significant in MP2. For survivors alone 15 the changes were stronger and more significant, indicating perhaps that mortality and recruitment of trees in the period was offsetting accumulation on survivors. The most robust results, statistically, came from the randomization testing for numbers of lianas per tree, confirming the decline between 1988 and 2018 in MP1 and MP2, in the proportion of trees with high numbers of lianas per tree, and the corresponding increase in the proportion of trees 20 with very few lianas. This form of testing, taking independent tree samples at the two census dates (i.e. samples having different survivors), also removes 'survival bias' due to trees with heavy liana loads often dying along with their loads and the liana counts becoming zero (Visser et al. 2018a, b). This bias has, inter alia, important consequences for understanding how liana loads, or 'burdens', affect tree population growth rates (Muller-Landau and Vissier 25 2019). Randomization also avoided the more essential problem of a lack of independence caused by survivors being counted in the samples of both census dates. As to be expected for a long-lived structural parasite, liana numbers increased with the age of the surviving host unless they were shed. Even so, trees die for many other reasons, and infestation by lianas is but one contributor.

5
The weak, or lacking, differences in mean numbers of lianas between censuses in the Dipterocarpaceae in both plots may indicate that for this family, summed across the size classes, lianas were being gained at a similar rate to them being lost without undue mortality of large and very large trees. As discussed later, the interaction between lianas and the very largest trees points to an important process operating in this dominant tree family. Most of the 10 large to very large trees in the Danum plots are dipterocarps (Newbery et al. 1992(Newbery et al. , 1996. For the smaller Euphorbiaceae, however, there was little difference in mean number of lianas in MP1 yet substantial decreases in MP2. The negative binomial regressions for the changes in numbers of lianas per tree, under randomization testing, furthermore highlighted the importance of the year x size interaction, which again differed between plots. The When considering rgr response, the causal direction is not unambiguous. On the one hand, high liana loads might have reduced rgr more than did low ones but, on the other, trees that were growing fast, may have discouraged the survival of lianas more than those growing 10 slowly. This might have been connected with branch shedding: faster growing trees would probably also shed branches faster than slower growing ones. If slower growing trees suffered from higher liana loads they would have been more likely to die than those faster growing ones with lower loads. When the potential for a tree to grow tall was high, the relative costs of branch shedding would be low, it is 'affordable', compared with the 15 converse. An important qualification to the analysis is that the rgr models were based on the growth of final survivors over the 30-yr period, and not on the rgr they attained up until the time they died on a shorter interval-wise basis using the tree census data (cf. Newbery and Ridsdale 2016).

Tree-liana interactions
Loss of lianas is closely connected with death of their hosts. Tree dynamics found using dates near to the two liana censuses can be more accurately derived by averaging over to the four periods between the main tree censuses of 1986of , 1996of , 2001of , 2007of and 2015of (Newbery et al. 1999of , 2011and unpubl. data). These further estimates reduce the length of interval bias inherent to mortality rate calculations (Sheil and May 1996). The first and last tree censuses were 2-3 years prior to the liana ones, and showed slightly lower differences in rates between plots compared with those based on the liana-censused trees. For all trees, MP1 had a 33-35% higher annualized mortality rate (ma) than MP2 across three size classes 30 -< 5 60, 60 -< 120 and ≥ 120 cm gbh (Appendix 2, Table S2). For medium-sized trees, this proportional difference between plots was similar for the Dipterocarpaceae, but it was > 60% for the Euphorbiaceae. Plots differed hardly at all for large and very large dipterocarps (1-8%) but were somewhat lower for large euphorbs (20%). (There were very few euphorbiaceous trees in the very large size class to allow estimates.) in respect to plot differences. Nevertheless, the apparent trends among the very largest trees in number of lianas with increasing gbh (Fig. 3, Appendix 2, Fig. S2a, b) is due to three liana counts of 20 and 18 (MP1), and 13 (MP2): one tree dying, two surviving. In both plots there 5 were only three trees with ≥ 7 lianas each.
Whilst tree mortality rates were higher in MP1 than MP2 for all three size classes, it was the disproportionately higher-laden trees in the very large size class that led to more lianas being lost in MP1 than in MP2 in that size class. The large and very large trees in MP1 most likely had higher liana loads than those in MP2 even before 1988. In MP1 they had then 10 more lianas to lose by 2018 than those in MP2 by 2018. Without precise data on large branch fall rates, the large to very-large trees must be assumed to have been shedding their lianas at similar rates in MP1 and MP2. Therefore, the difference between the plots, for a majority of the dipterocarps at least, was a result of tree mortality over and above a higher branch fall loss in very large trees compared with large and medium-sized ones.

15
For trees overall to show similar changes in distributions of numbers of lianas per tree in both plots requires a different process to have been operating in MP2 from that in MP1.
Here the main changes were in the medium-sized understorey trees, where tree mortality was much higher in MP2 than MP1. However, survivors in these size classes were also slowly accumulating lianas, may be as newly rooted ones or coming from fallen trees (Putz 1984).

20
Recruits, being younger and smaller than other trees, would have had little time to accumulate lianas, and hence their loads were relatively low by 2018. The different dynamics in the two plots resulted in similar frequency distributions of lianas per tree though. The overstorey dipterocarps seemingly have more the option of branch shedding to lose lianas, but the understorey euphorbs far less (Campbell and Newbery 1993). The euphorbs can for the most part only lose their lianas on the death of the host. An interesting feature is the higher proportion of trees with no lianas in MP2 than MP1, particularly for the dipterocarps.
If fast-growing medium-sized dipterocarps (Newbery et al. 1999), moving from under-to overstorey, were losing branches they might well have shed their lianas along with 5 them. Shedding lianas with branches is one obvious way to relieve the load and recover rgr for onward growth for an overstorey species, and thereby increase its survival under competition with not just the lianas but also neighbouring trees. As understorey species do not grow into the overstorey at all, these considerations are not so important. Understorey trees incidentally provide a means for lianas to step-up onto neighbouring overstorey trees. If 10 large trees remained infested in their canopies and yet have no lianas on their own stems at near-ground level, these lianas are presumably coming across from neighbouring mediumsized trees. This would cause a mismatch between stem counts and liana effects.
The influence of spatial autocorrelation on the fitted values of the independent variables' coefficients was often quite small, both when the dependent variable was number 15 of lianas per tree (regressed on tree size and elevation) and when it was tree response as growth or survival (regressed on number of lianas per tree and tree size). This indicates that spatial aggregation of numbers of lianas per tree was low and near to random, i.e. locally trees with many lianas and trees with few lianas were well mixed and not forming patches of similarly high or low liana-laden trees in the forest. This form of (spatial) aggregation is not   This long-term forest process was probably continuing to converge between 1988 and 2018, with plot-level tree basal area moving towards its maximum faster in, or more ahead of, MP2 than MP1. Mature closed forests approaching, or already at, the site-defined maximum carrying capacity (Newbery et al. 1992) might be expected to lessen their liana loads because 5 as the forest became denser in basal area terms, branches that fell would be less easily replaced on the large trees in the overstorey, and lianas would have increasingly less favorable light conditions in which to recruit and grow on the small and medium-sized trees.
Heavy liana cover in the understorey would further limit light reaching tree leaves, and hasten their decline under competition with less infested neighbours (Schnitzer and Bongers 10 2011). By contrast, any large branch fall also perturbs the canopy by causing small openings, and these let in more light to the understorey and hence promote tree and liana growth.
Species with high mean numbers of lianas per tree in 1988 had much larger decreases by 2018 than those with low numbers per tree (e.g. from 3.0 vs 1.0 to ~1.9 vs 0.9). This indicated a form of negative density dependence at the community level, perhaps by the more 15 heavily laden species disproportionally shedding more lianas or dying from them, or recruiting them more slowly. However, species' numbers of lianas per tree in 1988, in 2018, and the differences between censuses, were not correlated with mortality (%d) or growth rate of stems (rgr). Tree recruitment into the ≥ 30 cm gbh population might have been expected to have been highly correlated with growth rates of small trees (10 -< 50 cm gbh). This suggests Another reason why MP1 had more lianas on its very large trees than MP2 in 1988, 25 might have been topography (see Fig. 1 in Lingenfelder and Newbery 2009). MP1 is on the side of a long pronounced ridge with steep slopes and is relatively topographically distinct, whereas MP2, also with an overall similar difference in elevation of 35 -40 m over the 400 m plot length, is less steep and forms a spur to a higher backing ridge. These subtle differences in topography likely played a role in water availabilities during dry periods, the 5 top of MP2 being less prone to dryness than the ridge top in MP1 (N. Chappell, pers. comm.).
Past drought effects could ex hypothesis have been stronger in MP1 than MP2, making larger lighter areas in the understorey, and this promoting the establishment of more lianas in the former than latter, with lasting consequences.
Although the number of lianas decreased with elevation there was a similarly strong 10 relationship in both plots: so differences in elevation -at least in recent decades -were not as such a good explanation for differences in the different liana changes across plots. Very large trees were also quite evenly spread across the topographic gradient in both plots (D. M. Newbery, unpubl. data). The proportions of the very largest trees (gbh ≥ 240 cm) dying in the lower and upper halves of the plots (corresponding roughly to lower slopes and ridges) were 15 35 and 32% in MP1 and 46 and 50% in MP2 (χ 2 = 0.05, P = 0.82; χ 2 = 0.09, P = 0.77). So host size distribution was not an explanation either. This would then indicate that MP2 probably recovered faster than MP1 in successional terms and any very large trees that did have high loads in MP2 had already fallen before 1988, and in the 1988-2018 period a similar phase was being realized in MP1. Meanwhile MP2 was showing a next, maybe final, phase of 20 liana loss within the understorey. Indeed, if intensity of disturbance was less in MP2 than MP1, a faster recovery (less impeded by lianas) would have been expected more in the former than the latter. This represents a shifting mosaic in time and space, in theory reminiscent of the 'pattern and process' idea of Watt (1947), applied by Richards (1952Richards ( , 1996 to rain forests. The dynamics of lianas in such a late successional forest can be linked historically to liana activity soon after that distant strong disturbance set the long-term regrowth on its variable course.     App. 2. Table S1. Liana infestation on the eight other families (extension of Table 2). Table S2. Tree mortality rates by plot, period and main families.  Tables S5-S6. Interpretation of date-size interaction, and model predictions. Table S7. Regressions of number of lianas on date and size (model 4) for the two 10 main families (complementing Table 8).
App. 4. Table S1. Statistics, and model fitting, for the four compared models of number of lianas on date and size in Table 8. Table S2. Comparison of best fit models under randomization testing.
App. 5. Table S1. Liana infestation on the most abundant species, individually, at two dates.