Until the beginning of the 20th century, the American chestnut (Castanea dentata Borkh.) was a hallmark tree of eastern American hardwood forests, ranging from Ontario to Alabama and spanning from the Atlantic coast to Illinois (Russell, 1987; Collins et al., 2017). Throughout this range, C. dentata had crucial ecological, economic, and cultural importance–providing a valuable nut crop for wildlife (Diamond et al., 2000), rot resistant and durable timber for manufacturing (MacDonald et al., 1978), and properties enabling the ways of life of Native Americans and Appalachian communities (Steiner and Carlson, 2006). In 1904, an invasive chestnut blight, cased by the ascomycete fungus Cryphonectria parasitica (Murr.) Barr, was discovered on C. dentata trees (Rigling and Prospero, 2018). This fungus was unintentionally introduced to eastern American forests prior to 1904, probably on nursery stock from Japan (Milgroom et al., 1996; Dutech et al., 2012; Rigling and Prospero, 2018) and the blight spread rapidly, functionally extirpating C. dentata from the overstory in just 50 years (Paillet, 2002).

Since the loss of C. dentata from the forest overstory, considerable efforts have been made to introduce genes for blight resistance via introgression from the Asian species of Castanea into locally adapted populations of C. dentata throughout its native range. The American Chestnut Foundation (TACF) is piloting per-state chestnut backcross breeding programs to develop a hybrid American chestnut tree with C. dentata traits but blight resistance from other chestnut species, including C. crenata, C. henryi, and C. mollissima. Researchers at SUNY-ESF have independently been developing a transgenic method enhancing blight resistance in C. dentata (Steiner et al., 2017); recently, both TACF and SUNY ESF have converged these methods into one united approach (Westbrook et al., 2019). Both methodologies seem promising, but both require finding genetic material from surviving mature C. dentata trees. Field breeding methods are important, but it is also crucial to understand C. dentata habitat preferences to determine where to plant blight-resistant trees in the future.

Species distribution models (SDMs) are useful tools to predict probable areas of species presence as well as areas of suitable habitat for a given species and can contribute valuable knowledge on species extent across a landscape (Elith and Leathwick, 2009). By using SDMs over large areas, managers can efficiently isolate the most ideal locations for C. dentata habitat before using boots-on-the-ground approaches such as soil samples to choose the best sites for restoration. SDMs use layers of environmental data, species occurrence points, and species absence points to generate statistics and predictions of species distribution (Franklin, 2010). Environmental variables used for SDM are often layers describing the topography, land cover, climate, or soil attributes of the region that may impact suitable habitat. For example, historical accounts of C. dentata in Pennsylvania suggest that chestnut was typically found on sandy soils and ridge topography, which describe important environmental layers to include in SDMs (Nowacki and Abrams, 1992).

Species distribution modelings have been used to model habitat distribution for a variety of tree species in order to inform land management strategies in the face of climate change (Booth, 2018). Matthews et al. (2011) examined 134 tree species responses to climate change using SDMs and found that species life history characteristics played a role in range shifts. Previous research has also used SDMs to explore C. dentata distribution at various spatial scales and extents across its range. Full-range studies found that temperature, precipitation, and soil factors influenced C. dentata distribution (Barnes and Delborne, 2019; Noah et al., 2021). Finer-scale SDMs for individual states or sub-regions often identified soil and topographic variables as most influential. Fei et al. (2007) modeled habitat of American chestnut in Mammoth Cave National Park using ecological niche factor analysis and found that ridges and steeper slopes were strong predictors of chestnut habitat. Tulowiecki modeled the range of American chestnut trees in western New York using historical tree records and nine different SDM techniques [artificial neural networks (stocktickerANN), classification tree analysis (CTA), flexible discriminant analysis (FDA), generalized additive models (GAM), gradient boosting models (GBM), generalized linear models (stocktickerGLM), multiple adaptive regression splines (MARS), maximum entropy modeling (MaxEnt), and random forests (RF)], and found that soil pH and slope were important habitat predictors (2020). These multiple SDMs differ in their approach to modeling species distribution by whether they utilize statistical or machine-learning methods, whether they model linear and/or non-linear relationships between variables, and whether they can model interactions between variables, and comparison between approaches can enhance the reliability of results (Tulowiecki, 2020).

This study contributes to understanding of SDMs of C. dentata across the state of Pennsylvania, which is central to its former range, by identifying the strengths and weaknesses of different SDM approaches and enhancing model robustness through ensemble modeling. Results of this study can be used to find further genetic material for development of blight-resistant American chestnut trees and help identify locations for pilot restoration projects to focus their money and energy (Fei et al., 2012). Using an ensemble of species distribution models, following the methodology outlined in Tulowiecki (2020), this study aims to better understand the distribution of C. dentata across environmental ranges in the state of Pennsylvania, identify areas most suitable for reintroduction efforts, and determine how different modeling techniques perform in modeling spatial distribution of American chestnuts in Pennsylvania.

We used a dataset of mature C. dentata locations maintained by TACF to model species distribution across all of Pennsylvania. We utilized nine different SDM techniques in the “Biomod2” package from R statistical software utilizing the “ShinyBIOMOD” GUI (Thuiller et al., 2009, 2016; R Core Team, 2013). All environmental variables used for modeling were generated and processed using ArcMap 10.7.1 geospatial software (ESRI, 2011). Determining true absence points for the whole region was not feasible due to the scale of this study across the state of Pennsylvania so pseudo-absence points for C. dentata locations were selected randomly from the entire set of potential absences within the “Biomod2” modeling process (Phillips et al., 2009). These pseudo-absence points were used as absence locations for the SDMs that required these points.

Surviving C. dentata locations were acquired from TACF’s “dentataBase” of known and verified mature American chestnut tree locations (TACF, 2020). Most tree locations in this database have been collected by volunteers who are asked to send samples to TACF so that they can be verified by experts. These samples contain freshly cut twigs with mature leaves attached and the location of the tree. Presence points collected in this manner can contain uncertainty and sampling bias due to GPS inaccuracy and the greater likelihood of finding and marking trees near roads or trails where humans have access, but these points represent the most comprehensive dataset of verified surviving C. dentata locations and are thus the most useful input to SDMs despite their limitations. We filtered this database so that it would only include C. dentata records within Pennsylvania. Filtering the database resulted in 295 non-hybrid C. dentata points in Pennsylvania.

Ten environmental variables representing land cover, topography, and soil attributes were included in the species distribution models to predict C. dentata habitat suitability. Variables were identified to represent a range of habitat characteristics that could describe growing conditions for C. dentata based on prior knowledge of the species biology (Paillet, 2002; Collins et al., 2017) and other modeling research for chestnuts (Fei et al., 2007; Tulowiecki, 2020). These variables included canopy cover variables acquired from National Land Cover Database (NLCD) (Homer et al., 2012), soil composition collected from ISRIC soil grids (Poggio et al., 2021), a Euclidean distance to streams layer generated from an Environmental Resources Research Institute streams shapefile (Environmental Resources Research Institute, 1998), and seven digital elevation model-derived variables that represent a range of topographic and moisture conditions. The digital elevation model was acquired from a U.S. Geological Survey (2000). The digital elevation model-derived aspect layer was transformed using a Beers transformation to change a cyclic variable to a linear one for more accurate consideration in modeling (Beers et al., 1966). Other variables, such as soil pH, that have been considered important for chestnut habitat in prior studies were unavailable at the spatial scale and extent for this study and thus could not be included in our models. We acknowledge that our results may be impacted due to exclusion of such layers. All variables were resampled to a resolution of 238 m to match the resolution of the soil data. Table 1 shows all environmental variables used in the models, their value ranges, and their initial resolutions.

In preparation for SDM, environmental raster data layers were processed in ArcMap 10.7.1. All rasters were reclassified to have the same cell snaping, mask, and cell size (approximately 238 m²). We generated Elevation, Slope, Aspect, Curvature, Topographic convergence index, Topographic position index, and Topographic relative moisture index (Parker, 1982) layers from the Pennsylvania digital elevation model acquired from PASDA (Pennsylvania Spatial Data Access, 2022). We generated a Distance to streams layer by using the ArcMap Euclidian Distance tool on the streams shapefile acquired from PASDA. The species location points were generated from the TACF dentataBase and we recalculated their coordinates to match with the coordinate system of the environmental layers.

We used ShinyBIOMOD, a graphical interface for the R package “biomod2,” to streamline the SDM process. We defined a geographic region of the state of Pennsylvania and uploaded the species occurrence data and environmental data we previously generated to train and apply the models. As our dataset only contained occurrence records, we generated three pseudo-absence datasets of 290 randomly generated pseudo-absence points in order to run the SDMs. We ran nine different SDMs (summarized in Franklin, 2010) to model C. dentata distribution. These models were ANN, CTA, FDA, GAM, GBM, GLM, MARS, MaxEnt, and RF.

Presence and pseudo-absence points were randomly split into 80% training data and 20% validation data. We ran five replications of all models with each of the three pseudo-absence datasets (totaling 15 replications of each model) to determine C. dentata species distribution and model statistics. Each model replication recorded True Skill Statistic (TSS) and Area-Under-the-Curve (AUC) error evaluation metrics which are commonly used to evaluate model performance (Swets, 1988; Allouche et al., 2006; Fawcett, 2006). Finally, we built ensemble-models of all generated SDMs to evaluate C. dentata distribution as informed by all nine modeling approaches. We generated three different predictions of C. dentata distribution using ensemble modeling outputs. These outputs were predicted probability of presence, binary predictions of presence/absence, and predicted probability of presence by committee agreement of binary individual model outputs. Binary predictions were made by creating a threshold of predicted probability to maximize sensitivity and specificity of the model. This means that a threshold was applied to probability predictions to maximize the true positive rate (sensitivity) and true negative rate (specificity) across the entire range. Pixels with probability values below the threshold were categorized as absences and pixels with probability values above the threshold were categorized as presences.

Accuracy statistics for each species distribution model are summarized in Table 2. These values represent the mean values for AUC (area under the receiver operating characteristics curve) and TSS (True Skill Statistic) between all model replications. AUC values range from 0 to 1 and represents the probability that assigning a predicted suitability value at a random presence point is higher than assigning a predicted suitability value at a random absence point (Fawcett, 2006). Generally, AUC values between 0.6 and 0.7 are interpreted as “fair” models and AUC values between 0.7 and 0.8 are interpreted as “good” models (Swets, 1988; Fawcett, 2006). TSS ranges from −1 to 1, with 0 representing a model that performs no better than random guesses (Allouche et al., 2006). Aside from the random forest models which showed an AUC and TSS value of 1,000, the gradient boosting model and classification tree analysis had the highest AUC and TSS values. AUC and TSS values of 1,000 indicate perfect agreement or fit between errors of sensitivity and specificity (Allouche et al., 2006) and we suspect model overfitting occurred in the random forest model, skewing its results. Artificial neural networks and generalized linear models had the lowest AUC and TSS values.

We generated variable importance measurements for all environmental variables in each SDM technique to identify the most important predictors for C. dentata habitat (Table 3). Sand to clay ratio of the soil was the most frequent top predictor of C. dentata distribution, identified in eight out of nine models. These models identified a positive relationship between sand to clay ratio of the soil and probability of C. dentata presence (indicating higher probability of American chestnut presence in sandier soils). Canopy cover (identified as important in five of eight models) showed a positive relationship with probability of C. dentata presence (indicating higher probability of American chestnut presence in areas of denser canopies). Other important variables included topographic convergence index (three models with a negative relationship) and topographic position index (three models with a positive relationship). The negative relationship with topographic convergence index indicates lower probability of American chestnut presence in areas of higher water accumulation, while the positive relationship with topographic position index indicates a higher probability of American chestnut presence along peaks and ridgeline formations. The artificial neural network model differed the most from other models in identification of variable importance, indicating canopy cover, distance to streams, and elevation as the three most important variables for modeling C. dentata distribution.

We generated three maps representing predicted C. dentata species distribution across Pennsylvania through ensemble modeling of all SDMs. The first map (Figure 1A) displays the mean of SDM predictions with the influence of each SDM being weighted by the TSS of that model. This means that all models contributed to determining the suitability of each pixel, but models determined to be more accurate had more influence than models determined to be less accurate. This figure shows higher values of distribution probability in areas such as Pike County in the northeast of the state and along the Spine of Appalachia throughout the center of the state. High areas of C. dentata suitability can also be seen in Allegheny National Forest and along Lake Erie in the northwest of the state.

Figure 1B represents a model of binary presence/absence of C. dentata in Pennsylvania based on the combined ensemble model of SDMs. This model defined a suitability threshold in order to maximize the sensitivity and specificity of the prediction. All pixels with suitability values lower than the optimal threshold were defined as absences and all pixels with suitability values higher than the optimal threshold were defined as presences. We determined the optimal threshold to be a value of 0.58, which – when tested with the presence and pseudo-absence points – resulted in 74 true positives, 786 true negatives, 75 false positives, and 23 false negative. According to this binary presence/absence, approximately 16,088 square kilometers of Pennsylvania are considered suitable for C. dentata occupancy.

Finally, we generated a model of predicted spatial distribution of C. dentata in Pennsylvania determined by committee agreement of individual models (Figure 1C). This model was generated by first creating binary presence/absence models for each individual SDM technique, applying a threshold to maximize sensitivity and specificity. Each pixel of the final model is then classified based on how many of the individual SDM techniques classify it as a presence or absence. If 0 models classify it as presence, there is a high consensus of absence. If 1–3 models classify it as presence, there is low consensus of absence. If 4–6 models classify it as presence, there is a low consensus of presence. If 7–9 models classify it as presence, there is high consensus of presence.