Incorporating interspecific competition into species-distribution mapping by upward scaling of small-scale model projections to the landscape

There are a number of overarching questions and debate in the scientific community concerning the importance of biotic interactions in species distribution models at large spatial scales. In this paper, we present a framework for revising the potential distribution of tree species native to the Western Ecoregion of Nova Scotia, Canada, by integrating the long-term effects of interspecific competition into an existing abiotic-factor-based definition of potential species distribution (PSD). The PSD model is developed by combining spatially explicit data of individualistic species’ response to normalized incident photosynthetically active radiation, soil water content, and growing degree days. A revised PSD model adds biomass output simulated over a 100-year timeframe with a robust forest gap model and scaled up to the landscape using a forestland classification technique. To demonstrate the method, we applied the calculation to the natural range of 16 target tree species as found in 1,240 provincial forest-inventory plots. The revised PSD model, with the long-term effects of interspecific competition accounted for, predicted that eastern hemlock (Tsuga canadensis), American beech (Fagus grandifolia), white birch (Betula papyrifera), red oak (Quercus rubra), sugar maple (Acer saccharum), and trembling aspen (Populus tremuloides) would experience a significant decline in their original distribution compared with balsam fir (Abies balsamea), black spruce (Picea mariana), red spruce (Picea rubens), red maple (Acer rubrum L.), and yellow birch (Betula alleghaniensis). True model accuracy improved from 64.2% with original PSD evaluations to 81.7% with revised PSD. Kappa statistics slightly increased from 0.26 (fair) to 0.41 (moderate) for original and revised PSDs, respectively.


Introduction
Predicting the natural distribution of species across the landscapes has often focused on abiotic-centric species distribution models. These models often ignore biotic interactions and assume that such biotic factors at the scale of the forest stand or patch apparently averaged out, a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Study area
The Acadian Forest Region of eastern Canada [27] includes the three Canadian Maritime Provinces: Nova Scotia (NS; excluding the Cape Breton Highlands; [28][29]), Prince Edward Island, and all but the northwestern corner of New Brunswick. This investigation uses the Western Ecoregion of NS as primary study area. Ecoregions are ecological land classification units that delineate macroclimatic differences at a provincial scale [30]. The Western Ecoregion extends from Yarmouth to Windsor, including the Halifax peninsula. Geographically, the area is located between 43˚27 ' to 44˚56 ' North latitude and 64˚02 ' to 65˚47 ' West longitude, with a total land area of 16,904 km 2 , representing about 31% of the total provincial land base. Regional variation in climate is largely influenced by the area's proximity to the Bay of Fundy in the north and the Atlantic Ocean from the northeast-southwest. The region's climate is characterized by cold winters and warm springs and summers.
Regional forests are home to about 32 tree species and remain a diverse mix of both conifer and broadleaf species [28,31]. However, as a demonstration of the procedure, we apply the calculations to 16 target tree species from the 32 common species. These target species are sufficiently abundant to be modeled (>200 observations) and include eight conifer species (i.e.,

Data to generate maps of PSD
Data required to generate maps of PSD, include (i) a Digital Elevation Model (DEM) of the Western Ecoregion of Nova Scotia, (ii) precipitation surface, (iii) spatial variability of PAR and SWC generated from the Landscape Distribution of Soil moisture, Energy, and Temperature (LanDSET) model [9,10], and (iv) remote sensing-based calculations of growing degree days (GDD) [32]. Precipitation data were obtained from 25 Environment Canada climate stations in the Western Ecoregion of Nova Scotia [31]. The DEM was generated from 3-arc second resolution point-data (~70 m at 45˚N latitude) acquired from the NASA Shuttle Radar Topography Mission. Net incoming shortwave radiation served as input to the calculation of photosynthetically active radiation (PAR); estimate of net shortwave radiation was generated by the solar radiation-module of LanDSET. The soil water content (SWC) was generated by the soil water balance module in LanDSET; the algorithm includes a decreasing soil moisture availability function (including evapotranspiration, percolation, and surface runoff) and provisions for the accumulations (precipitation, water flow from upslope region).
A GDD map from thermal remote sensing data was developed for the area based on the standard definition of GDD [32], i.e., where T avg is the average daily temperature, T base is a base temperature threshold set at 5˚C, and i = 1. . . n, where 1 and n represent the start and end day of the growing season. Remote sensing data used in the development of an enhanced GDD-surface (at 30-m resolution,

Species-specific response function and PSD calculation
Curves describing the relationships between species occurrence probabilities and abiotic predictor variables of PAR, SWC, and GDD are based on generic functions scaling species response values between 0 and 1, where 0 represents highly unfavorable growing conditions and 1, optimal growing conditions. Species-specific responses to PAR, SWC, and GDD assumed the following forms, respectively: where c 1 is a scaling factor, c 2 is the slope of the light response curve, c p is the light compensation point or PAR at which the amount of carbon dioxide released in respiration equals the amount used in photosynthesis and the amount of oxygen released in photosynthesis equals the amount used in respiration [ Table 1].
[33], where SWC min , SWC max , and ѱ (optimal SWC; with SWC min < ѱ < SWC max ) denote species functional parameters that define the shape of the soil-water response curve.

½8
[34], where GDD min and GDD max are the minimum and maximum GDDs and represent the northern and southern limits of species tolerance for GDDs. The abiotic-centric PSD (i.e., PSD original ) was expressed as multiplicative interaction of species environmental response [9], i.e., The PSD values range between 0-1, where 0 represents unfavorable site conditions (and, thus, potentially low probability of species occurrence), and 1, superior site conditions (and potentially high probability of species occurrence).
Incorporating interspecific competition into original PSD. Cumulative effects of simulated interspecific competition over a 100-year simulation period with JABOWA-III provided the best possible estimate of species composition and aboveground productivity at the end of the simulation period. It is this information that we use to construct expressions of speciesspecific competitive ratings with Eq 10. Before scaling up JABOWA-III simulations of interspecies competition to the level of the landscape, we first classified the area into 12 forestland types (distinct combinations of tree species a particular site can carry) using surfaces of species distribution patterns generated with the abiotic-centric PSD model (i.e., eq. [9-10,35]). A classification of the area was needed to determine the actual range of forest growing conditions present and also reduce the number of JABOWA-III simulations. The methods employed in the classification and validations of the forestland types are more fully described in [24].
We initialized the JABOWA-III model with tree and environmental data from forest inventory plots located in each forestland type; all inventory plot data were in GIS format, courtesy of the Nova Scotia Department of Natural Resources (NS DNR). Individual tree growth, establishment, and mortality were simulated in JABOWA-III taking into account inter-tree competition, and growth-related physical variables such as sunlight, accumulated GDD, and soil water and nutrient content (http://dx.doi.org/10.5061/dryad.tf7f7). As in Eq (9) and LanDSET Table 1 [16,34], c 1 is a scaling factor, c 2 is the slope of the light response curve, and c p is the light compensation point c GDD min and GDD max are the minimum and maximum tolerance limits of species, based on a review of the scientific literature (e.g., [34]) and references from the internet-the same parameters are used in JABOWA-III d 0 and 1 represent the soil's permanent wilting point and field capacity, respectively [33] ѱ denotes a species functional parameter that defines the shape of the soil-water response curve.

Species
doi:10.1371/journal.pone.0171487.t001 simulations, mean monthly precipitation and mean monthly temperatures are based on data acquired from Environment Canada climate stations within the study area. Thus, all elements, but the soil conditions and tree-growth-related elements of JABOWA-III are addressed within Eq (9) (abiotic-only model) but at different temporal resolutions. Replicate simulations of the same forest patch (50 replicates) up to 100 years each yielded averaged tree production and composition trajectories. Competitive species ratings was defined by the total aboveground biomass of target species in cohorts consisting of other species, as a proportion of the maximum aboveground biomass of that particular species under optimum growing conditions on a plot. Values of species' competitive ratings were then rescaled to derive relative competitive ratings (P k 100 ) according to the following: where k and ' represent functional dependence on species k (whereby k = 1, 2, 3. . .. m) and forestland type ' (' = 1, 2, 3,. . .., n), respectively. Species-specific competitive ability by forestland type was then used to weigh original PSDs (by way of Eqs (9) and (10)) and give revised PSDs as follows:

Model accuracy assessment
The study area has a system of randomly placed forest-inventory plots. The network of plots provides a nearly continuous forest inventory dating back to the 1960s, incorporating field observations of species composition, growth, and mortality. Observed distributions of tree species in over 1,240 forest-inventory plots across the study sites were used to investigate the degree to which original and revised PSD values reflected actual species occurrence. All inventory plot data were also in GIS format, courtesy of NS DNR. For purposes of assessment, one assumption was that raster cells (i.e., 70 m x 70 m) representing the landscape were considered to have the biophysical attributes needed for species occurrence, where original and revised PSD values were greater than 0.25. The analysis summarizes relative frequency of each of the four possible outcomes as follows: (i) predicted vs. observed species occurrence (this outcome demonstrates positive agreement between predicted species occurrence and plot observation); (ii) predicted vs. observed species absence (this outcome demonstrates positive agreement between predicted species absence and plot observation); (iii) predicted occurrence vs. observed species absence (this outcome is indeterminate with respect to model predictions, as species absence from plots may be the result of other forest-forming factors not addressed in the current definition of PSD, e.g., as a result of species migration, disturbance, and forest conversion); and (iv) predicted absence vs. observed species occurrence (the outcome demonstrates potential inaccuracies in modeled biophysical factors and/or associated species environmental response). Analysis of model accuracy was based on two metrics, namely overall accuracy and kappa statistic [36,37]. The overall accuracy is the proportion of correctly predicted observations of species' presence and absence, whereas the kappa statistic corrects the overall accuracy of model predictions by the accuracy expected to occur by chance. The kappa statistics range from zero (very poor model accuracy) to one (perfect fit between predictions and observations).

Species competitive ratings by forestland type
Simulated values of relative competitive rating (P k 100 ) by species and forestland types are summarized in Table 2. During the 100-year successional period, species projected to possess the greatest competitive ability within a given forestland type have P k 100 -values = 1, whereas those expected to be eliminated from the community have P k 100 -values = 0. Final stages of stand development in all forestland types were dominated by relatively shade-tolerant, long-lived species, such as sugar maple, American beech, eastern hemlock, red spruce, and balsam fir, with significant components of eastern white pine, black spruce, and yellow birch. In the absence of large-scale catastrophic disturbance, as is assumed in our study, final species associations in most forestland types resemble a few of the old-growth forest types remaining on the NS mainland [29].

Original vs. revised PSD model
The original PSD model (i.e., Eq (9)) predicted four conifer (balsam fir, black spruce, eastern larch, and red pine; Fig 1) and two broadleaf (trembling aspen and yellow birch; Fig 2) species to possess the greatest potential to occur across the area. This predicted species occurrence is indicated by the abundance of high-quality sites (represented by yellow and brown colors; accounting for >70% of the land base) in their respective PSD maps.
Four conifer (red pine, red spruce, white spruce, eastern white pine; not shown) and two broadleaf (American beech and red maple; Fig 2) species were predicted to possess intermediate potential to occur in the landscape, as indicated by the abundance moderately high quality sites (i.e., green colors; accounting for >50% of the land base) in the corresponding PSD maps. The original PSD for one conifer (eastern hemlock; not shown) and four broadleaf (red maple, red oak, sugar maple, and white birch; not shown) species revealed that they may not fare as well as other species, because of a narrow response function to the prevailing abiotic factors used in the calculation of potential distribution. This is demonstrated by the abundance of low   (10)) lowered the potential distribution of all species. However, the extent of downgrading varies among species based on individual competitive success within forestland type. The revised PSD map suggested three conifer (balsam fir, black spruce (Fig 1), red spruce (not shown)) and two broadleaf (yellow birch (Fig 2), sugar maple (not shown)) species experienced nominal decline relative to their original distribution. Potential distribution for one conifer (eastern hemlock; Fig 1) and five broadleaf (trembling aspen and American beech, red maple (Fig 2), white birch, red oak (not shown)) species were projected to experience the largest reduction in PSD with the introduction of interspecific competition. For these species, most sites originally described high-to-moderate-quality diminished to low-quality sites, yielding over 86% increase in sites subjected to low levels of interspecific competition. Among all the 16 species studied, eastern larch (Fig 1) and trembling aspen (Fig 2) recorded the greatest reduction in high-quality sites to low-quality sites after adding a competition model to the original PSD model. However, pockets of moderate-quality sites (represented by green colors; Figs 1 and 2) existing across their respective revised PSDs indicate areas where eastern larch and aspen are expected to do well, competitively. Model accuracy assessments. Both the original and revised PSD models (Eqs (9) and (11), respectively) recorded a very high overall accuracy on average (Tables 3 and 4). Addition of interspecific competition improved overall model accuracy from 64.2 to 81.7%. The analysis based on the kappa statistic indicates that there is a greater degree of correspondence between actual species distribution in sample plots and revised PSD, than there is with original PSD (Tables 3 and 4). The kappa statistics increased from 0.26 (fair) with original PSDs to 0.41 (moderate) with revised PSDs.

Discussion
The original PSD model integrated modeled species-specific response to largely modeled biophysical variables of incident solar radiation (and PAR), SWC, and GDD; the model had been used extensively in eastern Canada to assess tree species habitat suitability [9][10]38], as well as classify the landscapes into forestland type [24]. However, the impacts of other forest-forming factors beyond the three biophysical variables used in the original definition of PSD were quite large, causing model inaccuracies to be as high as 73% in some areas [9]. We demonstrate a simple procedure for combining results from competition-based forest gap models to the current definition of PSD at the landscape level, without additional parameterization of the original model (i.e., Eq (9)). The critical question is whether scaling up output of gap models, in general, would be useful in explaining plant species distribution. Results show that accounting for competition effects in the original PSD model provides a better prediction of species distribution, as suggested before [39]. The decline in potential distribution for all tree species in the revised PSD model suggests that tree species cannot be expected to occupy all area of the climatic space that they can tolerate, due to the impact of other factors associated with forest development processes, including intraspecific and interspecific competition [3][4][5][6][7]. We found that mid-to late-successional tree species, such as balsam fir, black spruce, red spruce, yellow birch, and sugar maple, will be least affected by the incorporation of some measure of competition effect into the original PSD model. This effect is confirmed by the nominal decline in moderate-and high-quality sites in their respective revised PSD model. Many of these species are shade tolerant and are unlikely outcompeted for sunlight and thus may occur in many abiotically suitable sites [23,40].
Impacts of adding gap-model-based outputs of competition to tree species distribution were much greater when early successional and shade-intolerant species were considered. For example, eastern larch or trembling aspen, according to the original PSD, can occupy about 86% of the entire area as climatic space that they can tolerate. Incorporating their ability to compete for resources in cohorts consisting of other species reveals that they may not fare as well as initially predicted. As stand-alone species, they grow on a great variety of soil conditions and have the widest distribution of any native tree species in North America [41][42][43]. Although eastern larch and trembling aspen have potential to tolerate a wide range of environmental conditions, realistically, they are present in low abundance and contribute only 3.3% to the total tree volume in the area [44], contrary to what was predicted. Lower real abundance compared with predicted abundance can be attributed in part to low shade tolerance and Table 3. Results of an accuracy assessment between original PSD (Eq (9)) values and plot observations. Class limits for the assessment scale are based on Monserud and Leemans [28], namely <0.20 (poor), 0.20-0. 40  shorter life spans, rendering them less successful in association with other species over a 100-year period. Apart from low competitive ability, other factors that could cause their actual distribution to be lower include increased fire frequency, intensive harvesting, clearing for agriculture, insects, and diseases [41][42][43][44]. We did not use the disturbance option of the JABOWA-III model, other than canopy gapforming and self-thinning processes. However, there are some residual effects of disturbance captured in the abiotic surfaces, in particular in the expression of temperature and GDDs. Long-term GDD surface is based on 3 years of MODIS images of temperature that we subsequently calibrated with GDDs calculated at climate stations. In general, surface temperatures measured from space or at climate stations, for that matter, possess characteristics representative of the underlying surface. For example, air above forested surfaces, because of high evapotranspiration, tends to be cooler than the air above cutovers, grass surfaces, or airport runways during the day. These temperature differences would clearly become part of the overall expression of GDD, and in that sense would incorporate the effect of disturbance to some extent, introducing a level of noise in the prediction of PSD. Because SWC is based on a water balance calculation and surface temperature (used in the calculation of evapotranspiration, a component of the water balance), the impact of disturbance is also present. However, since SWC is mostly a function of landscape position and flow routing (redistribution of surface and shallow subsurface water; [9]), the effect of disturbance through surface temperature is expected to be considerably smaller. As JABOWA models forest succession based on initializations with tree data from forest inventory plots, the results from JABOWA could also indirectly account for Table 4. Results of an accuracy assessment between revised PSD (Eq (11)) values and plot observations. Class limits for the assessment scale are based on Monserud and Leemans [28], namely <0.20 (poor), 0.20-0. 40  the effects of disturbance. In this study, JABOWA-III takes these disturbances as initial conditions and projects these into the future without additional disturbance. Despite differences in timescales between Eq (9) and JABOWA, merging information from both models (by way of Eq (9)) does not lead to inconsistencies as they are used to describe a different aspect of PSD, i.e., spatial variation in species habitat (site) suitability vs. a onetime assessment of species aboveground biomass based on 100 years of cumulative effects of interspecies competition and changes in forest composition. However, the JABOWA model also includes a number of abiotic parameters that are not available in the original PSD model. Several site-specific factors addressed in JABOWA (e.g., soil texture, nutrient content, etc.) were not included in the development of the original PSD model at the landscape level. It is implicit that a model with more predictor variables is likely to have better performance [8]. However, over a 100-year simulation period with JABOWA, the role of these abiotic factors in changing species composition plays out across time, and what is left at the end of the simulation period is the relative competitive ability of the different species.
Furthermore, these additional factors were treated as random errors in the original PSD. By including these factors we could have potentially increased model accuracy; however, values for these variables at the landscape level are not readily available. On the other hand, if plot data could have been used to generate these same species' competitive ratings, naturally the expression of interspecific competition and their accumulation on model output would have been much more exact. Once the above challenges have been addressed, individual-based simulation models are likely to provide outstanding tools for overcoming past limitations and will provide the means to make reliable and robust predictions of the potential distribution of species at the scale of landscapes [6].

Conclusions
In this paper, we approach the problem of integrating biotic interactions, such as interspecific competition, into an abiotic-centric species distribution model by conducting an upward scaling of results generated with a forest-gap model (i.e., JABOWA-III) from shifting-gap mosaics to the landscape level, by way of Eqs (10) and (11). The results of the revised species distribution analysis indicate that, at the spatial scale of this study, the potential distribution of 16 common trees species in the area is largely a reflection of species' individualistic response to climatic factors and interspecific competition. Overall, the impact of incorporating some measure of interspecific competition into species' potential distribution across the landscape was relatively low for late-successional species with high shade tolerance that are dominant across the landscape (e.g., balsam fir, black spruce, red spruce, sugar maple, yellow birch) compared with early successional and shade-intolerant species (e.g., eastern hemlock, American beech, white birch, red oak, red maple, and trembling aspen).
The revised PSD model with interspecific competition (Eq (11)) includes a number of abiotic parameters that are not implicitly addressed in the original PSD model (Eq (9)). Differences between the abiotic-centric and the competition models may, in part, be due to the different input variables used. In future work, the original PSD model will be coupled with those abiotic variables absent in the model (e.g., soil texture, nutrient content).
for providing forest-inventory data in order that we may assess the accuracy of the predicted PSD surfaces against independent field-based information.