The complex, geologically young mountain region of Central Sulawesi consists of the West Sulawesi plutono-volcanic arc, mostly dominated by granite and other acid plutonic rocks and the Central Sulawesi metamorphic belt (Geological Research and Development Centre 1993, Hall and Wilson 2000, Hall 2013). The wide variety of soil types is due to variation in parent material, topography, elevation and climate; Acrisols, Cambisols, Ferralsols, Gleysols, Lixisols, Luvisols, Nitisols, and Regosols have been reported for LLNP (Dechert 2003, Häring et al. 2005, Culmsee et al. 2010b, Leitner 2010). According to climate data from 12 stations in and around LLNP (period 2002-2008; Kreilein, unpubl.), the climate of the study region is perhumid with most rainfall falling during April-May and November-December and a slightly drier period from July to October. At elevations below 1500 m in the inter-montane valleys, 1-3 drier months (<100 mm rainfall) do occur (Gunawan 2006). Humid north-westerly winds prevail from November to May and drier south-easterly winds during the rest of the year. ENSO-related droughts occur at irregular intervals (Gunawan 2006, Wündsch et al. 2014). While the populated valleys surrounding the park have a long history of human occupation and land use dating back at least 2.000 years (Henley 2005, Kirleis et al. 2011), most of LLNP itself has only been slightly affected by human activities like swidden agriculture, hunting and small-scale extraction of timber and non-timber forest products (NTFP) until the end of the 19th century. Since then, human impact has steadily increased and conflicts over resources and conservation efforts have appeared (Adiwibowo 2005, Weber 2006, Wündsch et al. 2014, Biagioni et al. 2015a).
Tree inventories
Criteria for study site selection were the presence of primary forest without major natural and human disturbance on level terrain (<10° inclination) without groundwater influence. The minimum distance between sites was 1 km except for sites N1 and N2 (c. 250 m distance). Five sites (S1050, S1400, S1850, S1900, S2400) had previously been surveyed in 2006/2007 (Culmsee and Pitopang 2009, Culmsee et al. 2010b, 2011) and were here censused a second time. Areas with any signs of logging or rattan-extraction were excluded. However, two NTFP, damar-resin from Agathis dammara (see Appendix 4.5 for full species names) and agarwood from Gyrinops sp., are collected by local inhabitants at elevations of c. 1700 – 2000 m throughout the study region and might have a slight, but unknown effect on forest structure and species composition.
We established 13 inventory plots (see Table 1, 0.24 ha size, 40 m × 60 m, without slope correction) and divided each plot into a 10 m x 10 m grid using plastic poles and string. Within the plot, we permanently tagged and pre-identified all trees with dbh (measured at 1.3 m) ≥ 10 cm, recorded their position and collected information relevant for species identification (habit, presence of sap, bark characters, etc.). Tree height was measured with an ultrasound clinometer (Vertex IV with T3-Transponder; Haglöf, Langsele, Sweden) and diameter at breast height (dbh) with a diameter tape. Epiphytes and lianas were removed from the trunk before measuring
dbh. The diameter of buttressed or stilt-rooted trees was measured above the buttresses (if these were taller than 1.3 m) and the measurement height recorded. For the measurement of multi-stemmed trees and strangler-trees, we employed the protocol described in Culmsee et al.
(2010b). The same procedure was followed in a 5 m × 5 m subplot in each of the 24 grid-cells of 10 m × 10 m (altogether 0.06 ha per plot) for small trees (dbh: 2 – 9.9 cm). The survey data of all plots will be made available through the repository ForestPlots.net (Lopez-Gonzalez et al. 2009).
Table 1 Characteristics of the 13 inventory plots in old-growth mountain rain forests in Lore Lindu National Park, Sulawesi, Indonesia. Temperature and rainfall data modelled for the time period 1950–2000 from WorldClim (2014). Soils according to IUSS Working Group WRB (2014). Data on bedrock material from Geological Research and Development Centre (1993).
Plot Elevational
1990 23.0 Sideralic Cambisol on acid plutonic rocks
S0850 Colline 850 Tokepangana 1°36.9' S
120°04.4' E
Level terrace 1990 23.0 Sideralic Cambisol on acid plutonic rocks
S1050
Sub-montane 1050 Pono Valley 1°29.7' S
120°03.4' E Level terrace 1900 21.0 Sideralic Cambisol on metamorphic rocks
montane 1950 Pantakleabae 1°42.0' S
120°09.0' E Gently
plateau 2080 14.8 Folic Histosol on acid
plutonic rocks
To characterise the availability of plant macro-nutrients and the chemical status of the soil, we dug four soil pits of ≥ 50 cm per plot. In each pit we took 100 cm3-samples from the mineral soil at 10 and 40 cm depth and from the lowermost organic layer (OF and OH horizons) to obtain volume-related nutrient contents. The samples were air-dried and transported to Göttingen for
nutrient analysis. The pH of the fresh soil was measured in the lab at Tadulako University, Palu, in a suspension of 10 g fresh soil in 1 M KCl. The concentrations of total C and N were measured by gas chromatography (vario EL, elementar, Hanau, Germany) from samples dried at 70°C (N) and corrected for moisture content.
Since no limestone occurs as bedrock in the study region, we assumed that all carbon was of organic origin. The total contents of Ca, K, Mg and Na, Al, Fe, and Mn in the organic layer material was measured by HNO3 digestion and subsequent ICP-OES analysis (Perkin Elmer Optima 5300 DV), and the concentrations of exchangeable Ca, Mg, K, Na, Al, Fe, and Mn in the mineral soil by BaCl2 extraction and subsequent element analysis in the percolate by ICP-OES. The observed pH change during the percolation process was used to calculate the concentration of hydrogen ions at the cation exchangers. Cation exchange capacity (CEC; in µmolc g-1) was obtained by adding the charge of K, Mg, Ca, Na, Al, Mn, Fe, and H measured in the percolate; we defined base saturation as the percentage of K, Mg, Ca, and Na ions (expressed as µmolc) in CEC. One sample from S1200 showed unrealistically high values for Ca and was excluded from the analysis. In the additional soil pits, we collect soil samples from their main detectable horizons and described the soils according to FAO (2006). We then complemented the soil descriptions with literature information (Dechert 2003, Häring et al. 2005, Culmsee et al. 2010b, Leitner 2010) to identify the soil type according to the WRB system (IUSS Working Group WRB 2014).
Collection of specimens.
In the plots that were surveyed for the first time, we generally collected herbarium specimens from each individual. From morphospecies that were easily recognized in the field (e.g. palms, Dracaena spp.), we collected at least one specimen per plot. In the plots that had previously been surveyed, we collected only trees that were either in flower or fruit, had not been determined to species level previously (Culmsee and Pitopang 2009, Culmsee et al. 2011), or were recorded for the first time in our survey. In addition, we collected specimens of flowering and/or fruiting trees growing outside all plots. We used a telescopic pruner for smaller trees;
for large individuals, we employed a slingshot (Big Shot, Sherrill Tree, Greensboro, USA) and hand chainsaw attached to strings. In a few cases, we used rope-climbing technique to access trees that were difficult to sample otherwise. Specimens were conserved in 70% denatured alcohol in the field until drying at the facilities of Tadulako University in Palu. In total, 2156 numbers were collected with at least three duplicates each, of which 16% were fertile, representing 51% of all morphospecies. Full sets of collections were deposited at CEB and GOET;
duplicates of fertile specimens have been sent to BO and will be sent to K and L (Thiers 2018).
A small number was sent to the specialists of single families at E, KEP, MO and STU. See Appendix 4.2 for details on species identification.
Data analyses
All calculations were performed using the software RStudio, Version 0.99.491 (RStudio Team 2015) based on R, Version 3.2.3 (R Core Team 2015). First, we investigated the correlation
between plot elevation and seven soil parameters with Pearson's r using the rcorr function of the Hmisc library (Harrell Jr. 2015). We did not use precipitation data due to a lack of reliable data with high resolution at different elevations.
The diversity measures observed species richness (0Dobs) and effective number of species (1Dobs), are expressions of so-called Hill numbers (Hill 1973):, the latter incorporating abundance without favouring common or rare species (Jost 2006). Both measures are sensitive to undersampling, especially in species-rich ecosystems like tropical forests. We therefore also calculated the rarefied/extrapolated values based on a base sample size (BSS) of twice the number of individuals of the plot with the smallest number of individuals (BSS = 210; see Chao et al. 2014) and used the resulting effective number of species per BSS (1D210) as primary measure of diversity. The same procedure was applied for genera and families. We used the iNEXT package for R (Hsieh et al. 2014) for the calculation of diversity measures and rarefaction-extrapolation curves. The primary diversity measure was then applied as response variable in linear regression models with elevation and seven soil parameters as explaining variables. First, we ran bivariate linear regressions and sorted the environmental variables by their explanatory power (r²). Then, we excluded all variables, which were significantly (p ≤ 0.05) correlated with the most explanatory variable, and repeated the process until only mutually uncorrelated variables remained. The resulting environmental variables were used in multiple linear regression models with the diversity measures as response variable. The models were then simplified by a backward selection procedure using F-tests until obtaining the minimum adequate models. For the relative number of endemic species per plot, we ran logistic regression models (LRM) using the function glm with binomial error structure and logit link function in R with plot elevation as explanatory variable.
We conducted a literature search for plot-based tree inventories in Malesia and recorded their locality, elevation, the number of individuals, the number of species and environmental information, where available. We then chose all plots from sources with reliable species identification (incl. deposition of voucher specimens in herbaria), a full list of the recorded species and their abundances and a plot size of 0.1-1.0 ha (n = 38, see Appendix 4.3) for comparison with our 13 plots. Here again, we used the iNEXT package for R (Hsieh et al. 2014) to produce rarefaction-extrapolation curves. As BSS, we also defined twice the number of individuals of the sampling unit with the smallest number of individuals (n = 90).
To search for patterns related to community composition, we also calculated the abundance-based index of Bray-Curtis dissimilarity (Bray and Curtis 1957) for all pairs of plots and used the resulting matrix as dependent variable for a multivariate analysis of variance (MANOVA) using the adonis function in the vegan library for R (Oksanen et al. 2016) with single environmental parameters (Appendix 4.4, Table 1) as explanatory variables and a Monte-Carlo permutation test with 500 permutations. We then sorted the environmental parameters according to their explanatory power (highest r²) and assessed the autocorrelation of these parameters using Pearson's r. All environmental parameters significantly related to elevation, the variable with
the highest explanatory power, were excluded. The practice was repeated with the most explanatory of the remaining variables. The resulting, mutually independent variables, elevation and pH of the organic layer (pHO), were both used for a second MANOVA, again with 500 permutations. Furthermore we used non-metric multidimensional scaling (NMDS, function metaMDS in the vegan package) to graphically display the influence of elevation and pHO on tree species composition.
Since the elevational gradient was most pronounced for both diversity and composition, we attempted to define distinct elevational zones based on species composition. We used the matrix of Bray-Curtis dissimilarities for all plot-pairs to conduct a cluster analysis using the hclust function of the vegan package for R (Oksanen et al. 2016) with the average linkage (UPGMA) method and plotted the results in a dendrogram. Subsequently, we performed an indicator species analysis (Dufrêne and Legendre 1997, De Cáceres et al. 2010) using the indicspecies package for R (De Cáceres and Jansen 2015) and tested the statistical significance of the associations using a permutation test (999 random permutations) to show all species having a significant (p < 0.05) association with any elevational zone.
To explore the affiliation of the tree species in the plots with main tree guilds, we pooled all trees with dbh ≥ 10 cm and small trees (dbh 2 – 9.9 cm) recorded in the 5 m × 5 m subplots (0.06 ha per plot; small trees were not recorded outside of these subplots) and assigned them to two guilds based on the observed maximum height of the species relative to forest canopy height using our measured tree heights and literature data. The canopy height of a plot was defined as the mean height of the 10% tallest individuals. All species reaching a maximum height of ≥ 2/3 of the respective canopy height in any of the plots were classified as ‘canopy species’.
Species not reaching this threshold height but which attained heights > 25 m according to the literature (e.g. van Steenis et al. 1948-2014), were also scored as canopy species. All remaining species were assigned to the guild of understorey species’ assuming that they possess special adaptations to the forest understorey environment and are more restricted in their height growth.
We then assigned the species in the two guilds to 35 orders according to APG IV (The Angiosperm Phylogeny Group 2016) , calculated the relative abundance of each order and guild per plot, and took the mean values for each elevational zone. For the 10 families with highest number of tree individuals, we ran logistic regression models (LRM) using the function glm with binomial error structure and logit link function using plot elevation and guild as explanatory variables and relative abundance of the order as response variable. We simplified the models stepwise by removing first non-significant (p ≥ 0.05) interaction effects and then variables until reaching minimum adequate models with only significant terms remaining or further simplification causing a significant (p < 0.05) increase in deviance as measured by an F-test.