Study site and period

The study was conducted at the Paracou Experimental site (5°18′N, 52°53′W) in a lowland terra firme rainforest of French Guiana dominated by the Caesalpiniaceae, Lecythidaceae, Chrysobalanaceae and Sapotaceae families. The site is located on private land and rented by the research institute (CIRAD). CIRAD is given full permission by the owner to conduct studies on the land, including the one presented in this manuscript. No specific permission was required for the activities conducted during the study (in situ measurements and collections). Field collections did not involve any endangered or protected species.

The climate of the area is equatorial with a dry season from mid-August to mid-November and a rainy season, often interrupted by a dry period in March. The mean annual rainfall and temperature are 3040 mm and 26°C. Soils are mostly shallow ferralitic and developed on schists and sandstones that form small elliptic hills.

The experimental setting of the site consists in fifteen 9-ha permanent sample plots (PSP), 9 of which were submitted to three logging treatments of increasing intensity (T1, T2, T3), while 6 remained control plots in natural conditions (T0). Each 9-ha plot contains a core zone of 6.25 ha surrounded by a 25 m wide buffer zone. In 1984, all trees with diameter at breast height (DBH) 10 cm were localized in the 250250 m core zone of each plot and botanically identified. Since then, all plots have been surveyed annually to evaluate individual growth, as well as recruitment and mortality. Logging operations were conducted during two years between 1986 and 1988 at the scale of the 9-ha plots. The treatments involved selective harvesting for timber (T1, T2, T3: 10 commercial trees.ha⁻¹ DBH or 60 cm), and fuelwood (T3: 30 non-commercial trees.ha⁻¹ with DBH between 40 and 50 cm), and additional thinning of non-commercial trees by poison-girdling (T2: 30 trees.ha⁻¹ DBH cm, T3: 20 trees.ha⁻¹ DBH cm; [28]).

In this study, we focus on the consequences of disturbance on population dynamics in terms of indirect effects. Direct effects of selective logging and poison-girdling include tree damage and destabilization by neighbours loss. These effects were materialized by the death of wounded trees up to four years after operations stopped. Here, we are interested in the response of species to the renewed availability of resources (space, light and soil) after disturbance. In order to eliminate the direct effects of the treatments, we thus calculated the population change rates and demographic parameters starting four years after the end of logging operations, over the period 1992–2003. We considered all trees cm DBH occurring in the six most heavily disturbed plots, i.e. in treatments T2 and T3.

Focal Species

Fifty-three species were selected for this study; they were well-documented in floras, unambiguously identified by field workers, and abundant enough to supply relevant data (see Appendix S1). They represented 28 families with a dominance of Caesalpiniaceae (13 species), Sapotaceae (4), Euphorbiaceae (4), Clusiaceae (4), Lecythidaceae (3) and Myristicaceae (3), which reflects the local composition of the forest. The studied species cover a large range of ecological strategies from small and medium-sized heavy shade-tolerant species to small light-demanding species and large light-demanding canopy species, including three pioneer species. In order to characterize these strategies, we estimated species growth potential in control conditions. We measured growth potential ( for relative diameter increase rate) as the 95^th percentile of the distribution of change rate in relative diameter increment

being the DBH at year , and taken here between 2000 and 2003 [27]. Note that, because there is no direct link between population change rates and species s, we chose to not include the latter in the presented bayesian network analysis. However, we reported their values in table S1 in Appendix S1 and performed the network analysis including for comparison purposes (Figure S1).

Demographic Parameters and Population Change Rate

We calculated the population change rate over the study period, , aswhere and are the number of living trees at the beginning (1992) and at the end (2003) of the study period of length . In order to take advantage of the yearly surveys, we calculated recruitment, defined as ingrowth over 10 cm DBH threshold (), and mortality () rates following [29] aswhere is the number of trees respectively recruited (i.e. reaching 10 cm DBH) or dying between years and , and is the number of living trees in year . We preferred this approach to estimate demographic parameters because estimations based on the raw comparison of an initial and a final survey [30] cannot account for trees both recruiting and dying during the study period. We calculated mortality rates excluding trees that had been wounded or poisoned during silvicultural operations in order to avoid direct effects of disturbance and capture intrinsic specific responses to disturbance. Note that although populations change in number as a result of the demographic processes of recruitment and mortality, the definitions of the corresponding rates do not lead to being a simple combination of the demographic parameters.

Functional Traits

Leaf traits were measured on leaves collected in three undisturbed control plots of Paracou following standardized protocols [31]. As far as possible, local measurements in control plots have been preferred to values from global trait databases in order to obtain locally accurate estimates of the interspecific trait variation [32]. Leaves were taken in fully sunlit positions, except for three species of small size (Iryanthera hostmannii, Oxandra asbeckii and Sandwithia guianensis). Ten leaves in total were collected on 3 different trees of each species. Analyses were performed on complete leaves, with petiole and rachis, except when rachises were too hard to be ground (Carapa procera, Jacaranda copaia, Schefflera decaphylla, Simaba cedron).

We measured the specific leaf area (SLA) as the one-sided area of fresh leaves divided by the oven-dried mass. Leaf nitrogen and phosphorus concentrations (LNC and LPC, the total amount of N and P per unit of leaf dry mass) were measured respectively by gas chromatography performed on a Carlo Erba CHN-OS elemental analyser, and by the molybdate-blue method, using ascorbic acid as a reductant. Leaf carbon isotope composition () was measured using a Finnigan DELTA S isotopic ratio-mass spectrometer.

As an estimate of species potential size, we considered the six control plots of the site in 2003 for the reference distribution of diameters and calculated the maximum diameter (D_m) as the 95^th percentile of the left-truncated distribution, that is all diameters above the threshold of 10% of the species absolute maximum, as recommended in [33]. We estimated seed size as the ellipsoidal volume calculated from measured dimensions:

and being respectively the length, the width and the depth of a seed. Seed mass was highly correlated to seed volume for a subset of species, but we preferred the latter trait because it was available for all studied species. Finally, wood density (WD), the ratio between mass and volume of a piece of wood measured at 12% moisture content was compiled from a local database and literature. Trait data for the focal species are presented in Table S2 in Appendix S1.

Bayesian Network Analysis

In graph theory, relationships within a set of studied variables can be represented by a Directed Acyclic Graph (DAG). A DAG consists in a set of nodes () and directed edges between them representing the variables and the dependencies among them respectively: for instance. In the Bayesian Network (BN) framework [34], [35], a statistical model of these dependencies is entirely defined by the DAG and an associated probability that measures the likelihood of the represented dependencies given the data ().

In a DAG, the parents of a node are the nodes that connect directly onto , noted . When conditioned on its parents, a particular node is by definition independent of all other variables that are not among the parents. Within each parents–child cluster, the dependency of the child node on its parents is thus described by a conditional probability distribution, , where is the variable associated with node . The set of distributions for all clusters in the DAG allows to define the likelihood of the DAG given the data

Specifically, the conditional distributions are modelled as Gaussian linear models [36]:(1)where and are regression coefficients and is the conditional variance. We thus log-transformed some variables in the analysis (, , and SV) to ensure that their distributions approached normality (Table S2 & S3 in Appendix S1).

Conditionnal dependence or independence among the studied variables can be directly read from a DAG. In the simple example , both nodes and depend on , so that and are dependent when conditioned on , i.e. when one considers as fixed [36]. In the example , node depends on , and depends on both (indirectly) and (directly), but and are independent conditionally on , that is and can vary independently when is fixed. Finally, in the example , node depends on both and , but and are independent because no edge exists between them. In the present study, the DAG has an inherent hierarchical structure given by the considered ecological model (Figure 1). In order to impose the sense of some of the arrows, relationships were specifically banned from the network [36]: we did not allow arrows directed from to demographic parameters or traits, neither from demographic parameters to traits, or from traits to .

We used a score-based approach to infer the structure of the DAG from the observed data: the algorithm goes through the space of possible networks, scores candidates based on the likelihood and finds the best scoring network [36]. We used a weakly informative prior for the joint probability, thus assuming weak a priori knowledge about the network structure. After a likely network was found, we performed an intensive heuristic search, which consists in perturbing the DAG a large number of times and performing multiple random restarts using the new network. This approach insures that a global maximum score has been reached.

We then performed complementary tests to ensure that the inferred multivariate relationships implied by the network were statistically valid. We first tested the overall structure of the network using Fisher's C test as described by [37]. We then tested all inferred relationships implied by the DAG structure ([37], [38], see Table S4 in Appendix S1) using partial correlation coefficientswhere measures the correlation between variables conditionally on variables in , and is the correlation between corresponding variables. We used Pearson's correlation coefficient to ensure consistency with the linearity hypothesis of the BN approach, as well as Spearman's correlation coefficient. All inferred conditional independencies were eventually verified. Inferred conditional dependencies present in the calibrated network that were not significant with regards to both coefficients () were removed from the analysis to ensure a conservative approach. P-values were corrected for multiple testing. All analyses were performed with R [39] using the deal [36] and ggm packages for BN and DAG analysis.

Population Change Rate and Demographic Parameters

After disturbance, recruitment critically controlled species demography, while mortality had revealed lower influence (Figure 2). Disturbance by logging induced instantaneous demographic imbalance in tree populations, because they were either directly eliminated by selective logging or poison-girdling, or indirectly killed by falling trees or following disturbance-induced stresses, implying that most of the final tree populations however recovered from reduced numbers [28]. This has resulted in increased recruitment [27], as well as increased growth after disturbance [40] when compared to populations in control plots. Regarding mortality, logging effects were more contrasted when compared to control conditions, with decreasing mortality rates observed in heliophilous species and increasing rates observed in shade-tolerant species, which indicated a negative effect of open canopy conditions on survival in the latter (see Table S1 in Appendix S1). These findings are consistent with previously documented disturbance effects on dynamics in logged-over stands [41]: dipterocarp forests showed higher recruitment rates in treated plots compared to control plots several years after mortality rates had returned to pre-disturbance levels. In forests of Venezuela, process-based simulations showed that ingrowth dominated over mortality after disturbance, leading to positive demographic imbalance [42]. We note that the time lag between logging operations and the study period, the length of the study period and the size threshold used to sample populations all determine how populations were sampled and may have consequences on the patterns highlighted here. Recruits come from germinated seeds or were present as advanced regeneration that has grown up to the threshold before the end of the study period. Because of their very low inherent growth rate, it is likely that recruitment in the more shade-tolerant species should come from advanced regeneration. A longer period of censuses would allow us to test the hypotheses of (i) a delayed and higher recruitment rates for shade-tolerant species together with (ii) a concomitant decrease in recruitment of light-demanding, and usually small-seeded species.

Drivers of Demographic Parameters

About half of the interspecific variability in the recruitment rate was best explained by three functional traits, including leaf phosphorus concentration, wood density and seed size (Figure 2). In disturbed stands of the site, we know that a positive relationship exists between growth and recruitment rates [27]. Thus the revealed relationships between functional traits such as leaf phosphorus concentration and recruitment rate could be mediated by species growth potential. Phosphorus is a major component of plant growth [43], as it is present in nucleotides implied in metabolic pathways (e.g. adenosine triphosphate, ATP). Under well-developed vegetation and on ageing soils, weathering and biological mobilization decrease phosphorus availability, which may in turn limit plant growth and development [44], [45]. Nitrogen is also known to limit tree growth [43], however N:P ratios observed here revealed high values, from 16.9 to 38.5 (26.2 5.4), which may indicate limitation by P, combined limitation by both P and N or even N limitation [45]. Recent evidence however confirmed that phosphorus rather than nitrogen controled the productivity of Amazonian forests [46]. Here, our network analysis evidenced the central role of phosphorus in the system, suggesting that species best at mobilizing available phosphorus also better recruited in response to disturbance. Unexplained variability in recruitment rates proceeds from unobserved causes that are sources of variation in juveniles survival and growth, it relates to stochasticity in the effects of abiotic and biotic factors on these processes.

Wood density negatively influenced recruitment rate, which is consistent with the well-known negative correlation between wood density and growth [15], [33], [47]–[49]. High wood density implies high stem construction costs, mechanical stability and resistance to cavitation [14], [50], [51], as well as low photosynthetic rates [52]. By contrast, low wood density is associated with efficient resource acquisition, which allows species to reform population numbers rapidly after disturbance. At individual scale, wood density is negatively linked to the tree ability to accelerate its growth during early life stages [53], that is light-wooded species have the best potential to accelerate their growth at intermediate sizes in suitable conditions, resulting in a hump-shaped growth trajectory. On the contrary, heavy-wooded species with constrained construction potential generally show slow [12], but relatively constant, growth patterns [53].

Seed volume negatively influenced recruitment directly, suggesting that regeneration processes were involved in species response to disturbance: high disturbance intensity tended to preferentially favor small-seeded species. However, as found elsewhere, we note that correlations between demographic parameters and seed volume were weak [5]. In undisturbed tropical forests, seed size was negatively correlated to mortality [5]. Although not causally linked to mortality of large trees, this correlation can be interpreted as a consequence of the adaptive trade-off between tolerance to low resource supply and growth [11], [23]. Large seeds often occur in species with attributes related to shade-tolerance, whereas pioneer and light-demanding species disperse overall smaller and more numerous seeds [54].

Previous studies found that mortality rates were negatively correlated with wood density in undisturbed tropical forests [5], . In our case, wood density was unrelated to mortality rate measured after disturbance. We suspect that logging-induced disturbance has blurred connections between the mortality rate and wood density. During the study period, light-demanding species, which are also light-wooded, survived better in disturbed than in control conditions (Figure 2) due to higher light availability [27], while disturbance had lower effects on the mortality rate of more shade-tolerant species. It also appears that correlations between wood density and demographic parameters might be stronger when considered over larger spatial scales (but see [6], [56]) leading to wider variability in wood density in relation with, for instance, variation in soil fertility [57], [58].

Species of high stature showed significantly lower mortality rates than smaller ones, consistently with results observed in undisturbed forests [27], [56]. We suspect that disturbance reinforced this relationship in logged-over plots: species of small stature are often shade-tolerant and, although some of them benefited from canopy opening [40], the majority probably suffered from sun exposure. In logged forests, gap creation increases lateral lighting and temperature along the vegetation profile. Such conditions have also likely impacted species differentially according to their ability to optimize water use. Interestingly, mortality also varied with enrichment in heavy carbon isotope, . In sunlit leaves, relates positively to water use efficiency (WUE; [59]), which measures the ratio between biomass gain through assimilation and water loss through evaporation. At Paracou forest site, water availability is the main climatic driver of tree growth [60]. It has been shown that tall tree species preferentially extract water from layers below 100 cm depth while shorter tree species show broader variations in the depth of water uptake [61]. Here, mortality decreased in species with increasing values of carbon isotope composition (higher WUE), but the correlation was indirect: the relationship appeared mediated by species stature and leaf phosphorus concentration. For a given maximum size, and for species with similar LPC, susceptibility to drought as indicated by carbon isotope composition, was not found to influence tree mortality.

Our network analysis evidenced the role of four functional traits as independent sources of interspecific variability. More specifically, leaf nitrogen concentration, wood density, maximum diameter and seed volume appeared as exogeneous variables, that is, variables not explained by other variables considered in the analysis. This finding constitutes an indirect validation of the leaf-height-seed (LHS) scheme which considers three major independent axes of variation in plant ecological strategies [9]. Among those four exogeneous variables, the proxy for the “leaf axis” of the LHS was leaf nitrogen concentration. Its relationships with the demographic parameters however were weak for recruitment and only marginal for mortality. Furthermore, SLA correlated better with LNC than LPC, LNC was more strongly correlated to LPC than SLA, and LPC appeared more strongly correlated to recruitment than LNC. Altogether, these relationships resulted in LNC acting as an independent axis of variation, LPC being a direct driver of recruitment and SLA having no direct influence on species response to disturbance. Wood density and maximum diameter characterized two different aspects of the “height” axis [12] in relation with resource allocation and conservation (WD) and resource acquisition (). Eventually, seed volume independently positioned species along the seed axis of plant ecological strategies.

Relationships among Functional Traits

We found close positive associations among three major traits involved in the so-called “leaf economics spectrum” (LES; [16]): SLA, LNC and LPC. These associations are indicative of coordination among the different processes leading to CO2 fixation by leaves, as SLA relates to light capture, LNC to carboxylation capacity, and LPC to energy transfer during photosynthesis. High values of these three traits lead to high rates of photosynthesis and respiration [16], [62].

The relationships between leaf carbon isotope composition, and the other traits of the LES are not straightforward to interpret. Relationships with leaf nutrients have not been found to be consistent across studies [63]–[65]. As an integrative proxy of WUE, is strongly driven by stomatal conductance and maximum photosynthetic rate. Here, leaf was found to be positively related to LPC (directly), but not to LNC, which may indicate a control of WUE by species photosynthetic activity rather than water flows. If one assumes low variance in stomatal conductance across species, these findings suggests that higher LPC implies higher because of lower concentration of CO in leaves due to higher photosynthetic rates [66]. Such a relationship was not observed with nitrogen, which also suggests limitation of photosynthesis by P rather than N.

A low water use efficiency (low leaf values) appears to be often related to a high SLA [63]–[65], [67], [68], but the lack of a relationship with SLA has already been found elsewhere [69], [70]. Clearly, a better understanding of how leaf relates to traits of the LES is required, which implies a thorough analysis of how carbon and water economies are coordinated within the leaf. We also found that leaf increased with species maximum diameter, a surrogate of maximum height. This pattern has been reported elsewhere both within [71], [72] and across species [20]. Gravity and xylem path length resistance induce a reduction in water potential at tree tops with increasing species height. This reduction is usually compensated for by decreased stomatal conductance leading to -enrichment and increasing [72], [73].