Reconciling the contribution of environmental and stochastic structuring of tropical forest diversity through the lens of imaging spectroscopy

Abstract Both niche and stochastic dispersal processes structure the extraordinary diversity of tropical plants, but determining their relative contributions has proven challenging. We address this question using airborne imaging spectroscopy to estimate canopy β‐diversity for an extensive region of a Bornean rainforest and challenge these data with models incorporating niches and dispersal. We show that remotely sensed and field‐derived estimates of pairwise dissimilarity in community composition are closely matched, proving the applicability of imaging spectroscopy to provide β‐diversity data for entire landscapes of over 1000 ha containing contrasting forest types. Our model reproduces the empirical data well and shows that the ecological processes maintaining tropical forest diversity are scale dependent. Patterns of β‐diversity are shaped by stochastic dispersal processes acting locally whilst environmental processes act over a wider range of scales.


survey.
We find that all three ways to calculate β-diversity from the field survey converge on similar estimates. Filtering under-story vegetation leads to estimates of β-diversity that are highly correlated to the estimates from the full survey (Mantel test,Mr = 0.99,p=0.001). The same is true for weighing trees' contribution to β-diversity by DBH (Mantel test,Mr = 0.99,p=0.001).
Supplemental Material, Figure S1: β-diversity calculated from the full field inventory survey (x axis) relates well to both β-diversity of canopy trees (left panel) and DBH-adjusted β-diversity estimates (right panel).

Hyperspectral sensor specification, calibration and data filtering
We used an airborne hyperspectral image of the Sepilok Forest Reserve obtained by a Specim AISA Fenix sensor. The sensor contains two photosensitive cells: one covering the visible and near infrared region and the other the shortwave infrared region to provide spectral range between 380 nm -2 500 nm. All wavelengths are captured by shared optics configuration and produce a single continuous image. VIS-NIR light is captured with spectral resolution of 3.5 nm and SWIR light with resolution of 12 nm. Further specification is available on the manufacturer website (https://www.specim.fi/products/aisafenix-hyperspectral-sensor/). The sensor was mounted onto a Dornier 228-201 aircraft operated by the NERC Airborne Research and Survey Facility (NERC-ARF).
NERC-ARF performed a wavelength calibration with known wavelengths emitted from Hg-Supplementary information, page 2 Ar, He, H, Kr, O, Ne and CO 2 lamps. A radiometric calibration was also performed by viewing an integrating sphere that meets the National Physical Laboratory standards. The integrating sphere provides known radiance at each wavelength that was used to produce a calibration curve applied to data collected with the sensor. Radiometric calibrations performed two months apart (December 2013 andFebruary 2014) demonstrated that the sensor is stable (average 0.15% change) with the exception of SWIR water absorption bands (2% change), a part of the spectrum not used in this study. The raw irradiance data (amount of energy reaching the sensor) was delivered to us by NERC-ARF alongside a data quality mask flagging bad pixels in the sensor (<1%) and synced but not applied flight information (i.e. data processing level 1). We applied the data quality mask and then georeferenced the image using the flight positional information, sensor field of view vectors and a high resolution digital elevation model of the tree canopies obtained by a ALS sensor during the same flight. As the field of view (FOV) of the hyperspectral sensor (18 • ) is wider than the FOV of the ALS sensor (12 • ), some of the image falls outside of the elevation model extent. We used the satellite-derived ASTER elevation model to patch the high resolution ALS model and provide a surface for the hyperspectral image to map onto. Pixels mapped onto coarse resolution satellite-derived surface models were later removed from our analysis as BRDF correction without good representation of the imaged surface is unreliable.
The mapped image contains artefacts due to scan and sun angles (BRDF effects) and atmospheric interference (aerosols, water vapour, etc.). We corrected for those artefacts using the MODTRAN radiative transfer model (Berk et al., 1999) as implemented in ATCOR v. 6.3.2.
The atmospheric correction standardised the raw irradiance values to reflectance values, making values from different flight lines comparable. We removed pixels illuminated over zenith angles over 50 • as BRDF correction might be unreliable. We also applied an NDVI filter (NDVI < 0.8) to remove pixels that do not exhibit a stereotypically plant signature, which might indicate they are coming from a non-forested patch or there has been a problem with obtaining or correcting the signal. Pixels with similar spectral signatures but different brightness might appear distinct on the image, hence we applied a brightness normalising function as described by Feilhauer et al. (2010). We also filtered spectral bands from each pixel that passed the quality filtering criteria above: the bands near the extremes of the sensor ranges are characterised by lower sig-Supplementary information, page 3 nal to noise ratio (wavelengths below 420nm and over 2400nm) and water absorption bands (1350 to 1480 nm and 1780 to 2032nm) are influenced by atmospheric water vapour. Cells centred on a particular wavelength are influenced by adjacent wavebands; we took advantage of this autocorrelation between adjacent bands to spectrally resample the pixels to preserve the signal while reducing the random fluctuations in the signal. We averaged three bands in the visible and the near infrared regions and two bands in the shortwave infrared region. Finally, individual flight lines were combined into a single mosaic by averaging overlapping pixels that passed the quality filtering outlined above, cropping the image to the boundary of the Sepilok forest reserve to exclude urbanised areas and palm oil plantations at its edges

Selection of number of clusters for dimensionality reduction
We used a modified version of the k-means clustering algorithm (Sculley, 2010) to group pixels by similarity. The k-means algorithm divides the observations into k number of groups, such that within-cluster sum of squares is minimised and between-cluster sum of squares is maximised. We tested a range of values for k and noticed that within cluster sum of squares decreases as the value of k increases, up to 250 clusters, when virtually all the variation in the data is captured by the clusters. The value of 250 for k was chosen for our analysis (Fig. S2).

Environmental maps of Sepilok
We used airborne laser scanning to map environmental variation over the Sepilok landscape.
The laser scanning sensor sends pulses of light down the forest canopy and times their return.
The first returns come from the forest canopy and enable us to compute a digital elevation model; some pulses reach all the way down to the forest floor and enables us to map the topography of the forest floor. Subtracting the elevation model from the terrain model provides us with a model of tree heights. We used this data to compute a number of metrics as described in the main text and visualised below (Fig. S3).
Supplementary information, page 4 Supplemental Material, Figure S2: Mean within cluster sum of squares decreases with increase of the number of clusters up to 250 clusters (black circle), after which little reduction is observed with splitting the clusters further.
Supplemental Material, Figure S3: Maps of ALS-derived environmental descriptors Supplementary information, page 5

Environmental variability between and within cells
We used generalised dissimilarity modelling to relate environmental to compositional turnover.
We found that the environmental variables that had the greatest influence on compositional turnover were the ones that vary more between cells than within cells. The figure below (Fig.   S4) illustrates within and between cell environmental variation. The cumulative distribution of mean plot environment is represented by a black line -steep slopes relate to little difference between cells and shallow slopes to more variation between plots. Within plot variation is represented by horizontal lines on the background of each panel that correspond to the interquartile range of each cell -short lines reflect consistent environment and long lines heterogeneous environment within the 1ha cell.
Supplemental Material, Figure S4: Environmental turnover (E) relates to compositional turnover as described in our GDM models (a-c). Cumulative distribution of mean 1ha cell environment is represented as a black line for the three ALS-derived environmental metrics that related to remotely-sensed β-diversity across the three forest types (d-l). Within-1ha cell environmental variability is represented by horizontal coloured lines, relating to the interquartile range of each cell.

Supplementary information, page 6
Neutral modelling approach Neutral simulations usually involve the following steps: (a) a random individual from a community is chosen to die and leave an empty space; (b) the replacement is chosen either to be a new species to the community (with some small probability ν, the per capita "speciation rate") or selected from propagules of existing individuals (with probability 1 − ν). Which propagule establishes in the space is determined using a dispersal kernel, representing the increased probability of nearby parents to provide offspring to the empty space. We used a two-parameter fat-tailed kernel (Clark et al., 1999) for elevated probabilities of long-distance dispersal events occurring. Typical forward-time neutral models require many iterations to reach a dynamic equilibrium and quickly become computationally prohibitive at large spatial scales. We use a backward-time coalescence approach to avoid these issues (Rosindell et al., 2008) producing equivalent results to forwards-time models with several orders of magnitude reduction in computation time. Coalescence allows for the full spatial structure of Borneo to be incorporated within the model -not just the Sepilok Forest Reserve -providing a more realistic, spatially accurate metacommunity. A single well-mixed habitat cell in our model was approximately 20m x 20m and contained a varying number of trees, depending on the density of the forest at that point. We used a combination of 50 dispersal and 7 speciation parameters for two sets of simulations, which resulted in 700 simulated virtual censuses over Sepilok. The range of parameters considered was beyond what can be expected biologically; the 7 speciation parameters were: 1 −7 , 5 −7 , 1 −6 , 5 −6 , 1 −5 , 5 −5 , 1 −4 ; We used latin-square sampling so that our 50 dispersal parameter combinations were evenly distributed across parameter space: measured in cell widths, σ ranged from 0.1 to 16 and τ ranged from 0.1 to 2. We used the laser scanning map of Sepilok to delineate the number of trees over the landscape using the approach of Coomes et al. (2017).
The tree densities outside Sepilok were assumed to be the average of the tree density of each forest type.
We developed two types of simulations to compare with the observed landscape patterns in Sepilok. The interactions between trees were treated as entirely neutral in the first set of simulations with the whole of Borneo considered as a single arena where the offspring position is only limited by dispersal. We refer to those simulations as "forest-type-naive". In the second Supplementary information, page 7 set of simulations, we overlay forest type structure onto the arena and refer to them as "foresttype-aware". We used published soil maps of Sabah (from the Sabah Forestry Department) and Borneo (from FAO/UN) to define the extent of the forest types outside of our immediate study region (Fig. S5). The forest-type-aware simulations still treat interactions within forest type as neutral, but mixing between forest types is restricted. The ability of species to recruit outside their preferred forest type was calculated from 12-hectares of inventory data available for each forest type. A species with over 50% of its stems in a particular forest type was considered to have a preference for that forest type, and their mean abundance in each of the other forest types compared to the preferred site was used as recruitment coefficients. Parts of Borneo covered by environments not present in Sepilok were assumed to be a mix of the two non-preferred forest types when calculating recruitment. Overall, the forest-type-aware simulations incorporate forest type as an environmental filter in an otherwise neutral model.
Allowing non zero density of all three niches in all environments permits dispersal between patches on a particular forest type that would otherwise be unreasonably restricted. In reality Sepilok is sufficiently isolated from primary forests that dispersal from outside its boundaries is no longer realistic, but its current assembly may well have been influenced by historic large forest extents that existed before the establishment of intensive tropical agriculture.
Note that we made no attempt to incorporate finer niche structure in our models. Dispersal and speciation parameters are key to the structure of diversity in neutral simulations. However, processes that operate below the 1-ha scale can manifest as apparently random fluctuations in β-diversity that are incorporated into this parameter, so we make no attempt to interpret it in a biological context, instead focusing on spatially autocorrelated diversity and the dispersal parameters that can reconstruct it. Patterns of spatial autocorrelation in biological diversity can be caused by limited seed dispersal (D), which results in close-together cells having more similar species composition than those further apart. It is also affected by environmental factors that may themselves be spatially autocorrelated (e.g. forest types occupy large patches with similar elevation or terrain ruggedness) so close-together cells are more likely to have the same forest type and environments (E×D). In order to select the parameter sets that best fit the observed pattern, we compared GDM outputs between the modelled communities and the remotely-sensed β-diversity patterns. The neutral models reconstruct a spatially-explicit com-Supplementary information, page 8 munity that matches the boundaries of the Sepilok Forest Reserve. We can, therefore, analyse the virtual census using the same GDM approach as we did with the remote sensing data. GDMs output a set of splines that describe the contribution of the predictor variables to the measured dissimilarity; we compare those splines and select the models that most closely reconstruct the GDM models from the airborne survey. We compared the forest-type-aware models with the full extent of the spatially autocorrelated pattern (D + D×E) of each forest type. We repeated the model selection, but this time requiring the models to only reconstruct the spatial structure that does not covary with the environmental component of the models (D component only).
Supplementary information, page 9