Climate Change and Local Host Availability Drive the Northern Range Boundary in the Rapid Expansion of a Specialist Insect Herbivore, Papilio cresphontes

Wilson, J. Keaton; Casajus, Nicolas; Hutchinson, Rebecca A.; McFarland, Kent P.; Kerr, Jeremy T.; Berteaux, Dominique; Larrivée, Maxim; Prudic, Kathleen L.

doi:10.3389/fevo.2021.579230

ORIGINAL RESEARCH article

Front. Ecol. Evol., 02 March 2021

Sec. Conservation and Restoration Ecology

Volume 9 - 2021 | https://doi.org/10.3389/fevo.2021.579230

This article is part of the Research Topic Benchmarking Biodiversity in an Era of Rapid Change View all 9 articles

Climate Change and Local Host Availability Drive the Northern Range Boundary in the Rapid Expansion of a Specialist Insect Herbivore, Papilio cresphontes

$\r\nJ. Keaton Wilson*$ J. Keaton Wilson^1*

Nicolas Casajus²

Rebecca A. Hutchinson^3,4

Kathleen L. Prudic^1*

¹School of Natural Resources and the Environment, University of Arizona, Tucson, AZ, United States
²FRB-CESAB, Montpellier, France
³Electrical Engineering and Computer Science, Oregon State University, Corvallis, OR, United States
⁴Department of Fisheries & Wildlife, Oregon State University, Corvallis, OR, United States
⁵Vermont Center for Ecostudies, White River Junction, VT, United States
⁶Department of Biology, University of Ottawa, Ottawa, ON, Canada
⁷Canada Research Chair on Northern Biodiversity, Centre for Northern Studies and Quebec Center for Biodiversity Science, Université du Québec á Rimouski, Rimouski, QC, Canada
⁸Insectarium, Montreal Space for Life, Montreal, QC, Canada

Species distributions, abundance, and interactions have always been influenced by human activity and are currently experiencing rapid change. Biodiversity benchmark surveys traditionally require intense human labor inputs to find, identify, and record organisms limiting the rate and impact of scientific enquiry and discovery. Recent emergence and advancement of monitoring technologies have improved biodiversity data collection to a scale and scope previously unimaginable. Community science web platforms, smartphone applications, and technology assisted identification have expedited the speed and enhanced the volume of observational data all while providing open access to these data worldwide. How to integrate and leverage the data into valuable information on how species are changing in space and time requires new best practices in computational and analytical approaches. Here we integrate data from three community science repositories to explore how a specialist herbivore distribution changes in relation to host plant distributions and other environmental factors. We generate a series of temporally explicit species distribution models to generate range predictions for a specialist insect herbivore (Papilio cresphontes) and three predominant host-plant species. We find that this insect species has experienced rapid northern range expansion, likely due to a combination of the range of its larval host plants and climate changes in winter. This case study shows rapid data collection through large scale community science endeavors can be leveraged through thoughtful data integration and transparent analytic pipelines to inform how environmental change impacts where species are and their interactions for a more cost effective method of biodiversity benchmarking.

Introduction

Biodiversity benchmarking is fundamental to both basic and applied ecological research offering insights into the biological processes shaping species and their interactions. Benchmarking is a labor intensive endeavor, often limited by participation and training. Recent advances in sensing technology and communication have led to a diverse and plentiful data landscape coordinating and improving biodiversity community science efforts at scale so that they can be used in meaningful ways for benchmarking efforts (e.g., Sullivan et al., 2009; Prudic et al., 2017). Observational web platforms and smartphone applications, automated camera arrays, and machine learning-assisted identifications have also changed how biodiversity data is collected, processed, and verified (e.g., Sullivan et al., 2009; Prudic et al., 2017) although challenges remain (Bonney et al., 2009). These technologies have expedited the rate of understanding and changed the research focus to exciting new areas where an informatics toolkit is now a necessity (Feng et al., 2020). One new aspect of benchmarking biodiversity is to evaluate where species are and which species they co-occur with, or species distributions and their changing interactions (e.g., Bueno de Mesquita et al., 2016; Palacio and Girini, 2018).

Species distributions are known to be greatly influenced by climate (Brown et al., 2016). Climate-related range shifts have been and are continuing to be documented globally across taxa and systems: terrestrial (Parmesan and Yohe, 2003), marine (Poloczanska et al., 2013), and aquatic (Rahel and Olden, 2008). With current changes in global climate, species range shifts (Parmesan et al., 1999) and extensions in both altitude and latitude are being observed (Roth et al., 2014; Kerr et al., 2015). While many studies have examined the ongoing changes in climate and their effects on biodiversity and species ranges, most consider only abiotic factors in their analyses, missing the potential importance of local interspecific interactions once a species moves into a novel environment beyond its previous range (Blois et al., 2013; HilleRisLambers et al., 2013; Wisz et al., 2013).

Several interspecific interactions are known to play important roles in shaping range boundaries including competition (Connell, 1961; Huey et al., 2009; Stanton-Geddes et al., 2012), mutualism (Chalcoff et al., 2012; Moeller et al., 2012), facilitation (Bader et al., 2007; Stueve et al., 2011; Ettinger and HilleRisLambers, 2017) and natural enemies (Freeman et al., 2003; Speed et al., 2010). When a species extends into a new local environment, there are a few main scenarios it can encounter (Holt, 2003; Urban et al., 2007; Sexton et al., 2009): (1) ecological conditions are similar enough to previous conditions that there is little immediate effect on fitness and population growth rate, (2) the new local environment may possess biotic or abiotic conditions that differ from the original local environment and can accelerate (e.g., competitive or predatory release; or (3) decelerate (e.g., nutrient or nesting limitation) range expansion.

For insect herbivores, climate change can influence abundance and distribution through direct mechanisms (physiological impacts on growth, development and reproduction that impact fitness) and indirect mechanisms (impacting biotic factors such as host plant quality or predator abundance) (Bale et al., 2002; Deutsch et al., 2008; Kingsolver et al., 2011; Robinson et al., 2017). How and when climate change will affect herbivorous insect dynamics has received considerable attention generating a diversity of observed responses, especially in the pest management literature (Porter et al., 1991; Cannon, 1998; Harrington et al., 2001; Altieri et al., 2015; Castex et al., 2018). Some species are expanding in ranges and abundance (Battisti et al., 2005; Robinet and Roques, 2010; Robinson et al., 2017) while others are retracting and decreasing in numbers (Robinet and Roques, 2010; Zvereva et al., 2016; Sánchez-Bayo and Wyckhuys, 2019). Host plant abundance and distribution play a key role in generating these patterns as herbivorous insects are often limited by larval food resources (Dempster and Pollard, 1981; Pearson and Knisley, 1985; Ylioja et al., 1999). Exactly how host-availability translates into patterns of distribution, abundance, and range shifts for insect herbivores is still contentious and particularly complex when combined with direct effects on physiology (Louthan et al., 2015; Lany et al., 2018). Our understanding of the determinants regulating species distributions are becoming more nuanced as we begin to incorporate information on species’ dispersal capacity, population abundance trends, and climatic variables into our models (Elith and Leathwick, 2009).

In this study, we investigate the role of host availability and climatic variables on the range expansion of the specialist giant swallowtail butterfly (Papilionidae: Papilio cresphontes) in northeast North America over the last 60 years (1959-2018), with an emphasis on the perceived accelerated expansion of the last 18 years. We combine evidence from raw occurrence data with a series of species distribution models for P. cresphontes and associated host plants to evaluate the rate and direction of range changes in relation to both abiotic and biotic factors. While other studies have incorporated biotic variables as model inputs (Bueno de Mesquita et al., 2016; Palacio and Girini, 2018), our approach was to model the distribution of the insect herbivore and host plants separately and using these independent models to make post hoc inferences and comparisons of ranges. Because both this insect and its primary larval host plants (the common prickly ash [Rutaceae: Zanthoxylum americanum], southern prickly ash [Rutaceae: Zanthoxylum clava-herculis] and common hop tree [Rutaceae: Ptelea trifoliata]) are conspicuous, they are often reported in systematic biological surveys and museum collections. In this study, we bring together a combination of museum collection, survey, and citizen science data to understand how host plant availability, climate changes, and butterfly abundance are influencing the rapid expansion of an herbivorous insect as a case study. This study is one of few to demonstrate the interplay of both climate change and biotic interactions in shaping range limits while focusing on the ecologically important role of herbivores.

Materials and Methods

Study Region and Time Interval

We focused on eastern North America (study area bounded by −94° and −65° longitude and 25° and 55° latitude) where Papilio cresphontes has been reported to be expanding rapidly (Finkbeiner et al., 2011; Breed et al., 2012) and data are readily available for both P. cresphontes and larval host plants, (Zanthoxylum americanum, Zanthoxylum clava-herculis and Ptelea trifoliata). Though records of P. cresphontes exist further west than −94°, we set this cutoff to minimize complications of misidentification and complex species boundaries with its congener P. rumiko. We categorized and compared two time periods: T1 (1959-1999) representing the period prior to the beginning of the rapid range expansion and T2 (2000-2018) as the period when the rapid range expansion to the north began. This cutoff point was determined from raw occurrence data (Figure 1).

FIGURE 1

Figure 1. Evidence of northward range shift of P. cresphontes from raw occurrence data. (A) The maximum latitude of a recorded occurrence of P. cresphontes by year. Larger circles indicate years with higher numbers of occurrence records (high numbers in more recent years are due to increased citizen-science activity). The dashed lines represent the breakpoint between T1 and T2. Three years with extremely low maximum latitudes (<35) were omitted for clarity. The blue line and gray bar represent the loess smoothing curve and 95% confidence interval. (B) A ridge-plot of kernel density estimates of occurrences for P. cresphontes. Vertical dashed lines represent latitudes of major cities within the range. Years with <5 occurrence records were removed from plotting.

Data Sources

Butterfly and Host Plant Data

Papilio cresphontes (Papilionidae) is a sub-tropical butterfly widely distributed across North America. P. cresphontes and host plant occurrence data were obtained from a variety of sources: iNaturalist¹, n = 3,007, Global Biodiversity Information Facility (GBIF²), n = 14,181, the Maine Butterfly Atlas³, n = 11, the Maritime Canada Butterfly Atlas⁴, n = 6, Massachusetts Butterfly Club, n = 512, Butterflies and Moths of North America⁵, n = 1,188, and eButterfly⁶, n = 3,083. Data from iNaturalist and GBIF were downloaded using the spocc package for R (Chamberlain et al., 2016). We filtered iNaturalist data to include only research-grade records before combining with other data sets. Combined data were filtered for time frame, duplicates, and study area extent (see below) before further analysis and model building. In total, we used 8,051 occurrence records for P. cresphontes and 2,697 occurrence records (combined) for all three host plant species.

Environmental Data

We used the TerraClimate data set (Abatzoglou et al., 2018), a 4 km × 4 km resolution gridded set of monthly climatological data from 1958 to 2017 (at the time of writing this) to generate environmental predictor variables for modeling. We calculated a set of yearly summaries of 19 bioclimatic variables (Fick and Hijmans, 2017), frequently used in species distribution modeling, using the dismo package in R (Hijmans et al., 2017) for each year in each time period (T1 and T2) and then averaged these summaries across each time period to provide temporally appropriate climate summary for each set of models. We included all 19 bioclimatic variables as predictors for modeling.

Species Distribution Models

Distributions of P. cresphontes and host plants were estimated using MaxEnt 3.4.0, a machine learning algorithm based on the principle of maximum entropy (Phillips and Dudík, 2008; Elith et al., 2011; Phillips et al., 2017). MaxEnt is a presence-background method, which is considered to perform well when modeling climatic niches across a variety of sample sizes (Wisz et al., 2008). We used the ENMeval package for model-building, testing, and tuning (Muscarella et al., 2014), ultimately building 8 total models (P. cresphontes and three host species for each time period).

We used a combination of geographically structured and regular k-fold cross validation for model testing and tuning. We generated 10,000 random background points per species-time period combination (within the geographic extent outlined by the occurrences across both time periods – a rectangle defined by the minimum and maximum latitude and longitudes of occurence points) per model and used the blockCV package (Valavi et al., 2019) to divide our study area into 400 km × 400 km blocks. Blocks were randomly assigned to folds 1-5 over 250 iterations to determine a block design that maximized evenness of occurence and background points spread across all folds. This procedure was repeated for every model (8 times in total). Occurrence and background points from folds 1-4 were used as training data for MaxEnt cross-validation and tuning, while fold 5 was reserved as a set of out-of-sample test data for final model evaluation. Throughout the manuscript, we refer to these data as test data. We used another set of random fivefold cross validation within the training data to tune model parameters (within the ENMeval package). Throughout the manuscript, we refer to these data as validation data. We tested linear, quadratic and hinge features (and all combinations) as well as a set of regularization multipliers (0.5-4 in 0.5-step increments). We examined models using a range of evaluation metrics (Supplementary Figures 1–8), but eventually chose the model with the highest area under the receiver operating characteristic curve (AUC) on validation data. All evaluation metrics were reported for the separate set of spatially explicit test data generated by blockCV (Table 1). AUC values typically range between 0.5 and 1, and can be used for relative comparisons between models with the same data (with higher values closer to 1 indicating models with better predictive capacity (Lobo et al., 2008). Once the optimal parameters for a given species and time-frame were determined, we built full models using all available occurrence data to generate predictions for subsequent visualizations and analyses. We mapped the “cloglog” MaxEnt output, which can be interpreted as probability of occurrence under the assumption that the species presence or absence at nearby sites are independent (Phillips and Dudík, 2008; Elith et al., 2011; Phillips et al., 2017). Importance of predictors was assessed using the permutation contribution metrics generated when building full models. These metrics are built as MaxEnt steps through modifications of coefficients for single features. For each variable, values are randomly permuted on training data and a model is reevaluated on the permuted data. Then, the resulting drop in AUC scores are tracked and normalized to percentages (Phillips et al., 2006). Thresholds for binary presence-absence maps and presence distributions were generated using the maximum test specificity plus sensitivity (Liu et al., 2005). For all models, we used species-specific (but not time-specific) geographic extents during model building and tuning, as well as making predictions for graphical outputs. Kernel density plots are used to show latitudinal distributions of model predictions and northern range limits.

TABLE 1

Table 1. Model parameter set and evaluation metrics on geographically structured test data.

MaxEnt has become a popular modeling resource because of its predictive power, ease of use, and a well-detailed literature to get researchers started (Phillips and Dudík, 2008; Elith et al., 2011; Phillips et al., 2017). However, this framework has also received criticism, with researchers advocating for more explicit examinations of tuning parameters, evaluation metrics, and the incorporation of tools to deal with sampling bias (Radosavljevic and Anderson, 2014). Recent software additions have addressed some of these challenges, and opened up the “black-box” of MaxEnt (Phillips et al., 2017), though issues remain, particularly in the transparency of researchers’ hyperparameter tuning and evaluation (Morales et al., 2017). To this end, we implemented recently developed tools (ENMeval and blockCV packages in R; (Muscarella et al., 2014; Valavi et al., 2019) to explicitly outline tuning (Supplementary Figures 1–8), and to incorporate a spatially independent evaluation design to minimize overfitting (along with the built-in regularization in MaxEnt).

Northern Range Limits

We calculated the distance between the northern limit modeled for P. cresphontes for T1 and T2 using a longitude class approach (Leroux et al., 2013). For each 4-km longitude class (i.e., each “column” of 4 km of longitude across the entire study area), we determined the latitude of the northernmost grid cell where the species was predicted to be present during T1 and T2. We selected the latitude-pairs (pairs of data for a single latitude at T1 and T2) for which we had grid cells with occurrence for P. cresphontes in both time periods for each longitude class and tested whether the average northern limit distribution of P. cresphontes differed between T1 and T2, using a paired t-test. We used similar methods to determine differences between northern range limits of P. cresphontes and Z. americanum for both time periods.

Results

Evidence of Northward Range Shift of P. cresphontes From Raw Occurrence Data

Patterns of occurrence (as opposed to the predictive outputs from species distribution models) indicate a strong trend of a rapid and recent northward range expansion in P. cresphontes since the earliest recorded records of the species in our dataset (1959). The butterfly’s highest recorded latitude in a given year has increased dramatically since 2000 (Figure 1A), and the predicted suitability has shifted from low to high in many cities close to the current northern edge of the range (Figure 1B).

Predictive Accuracy of Species Distribution Models

Maxent models with optimal complexity settings were chosen via hyperparameter tuning, and a variety of evaluation metrics were calculated (Supplementary Figures 1–8), but ultimately the feature classes and regularization multiplier of the model with the highest average validation AUC was used for each species-time period pair. Once the final parameter set was chosen, models were evaluated on spatially explicit out-of-sample test data created by blockCV. Overall, models had high predictive accuracy on test data, with AUC scores ranging from 0.871 to 0.957 (Table 1). Generally, models were complex and incorporated combinations of feature classes paired with regularization multipliers (Table 1). Final models were generated using the parameter set (feature classes and regularization multiplier) described above, but built with the full set of data (training + test) to generate predictive maps (Figures 2, 3) and distributions (Figures 4, 5).

FIGURE 2

Figure 2. MaxEnt model predictions for the giant swallowtail butterfly, P. cresphontes for t₁ (1959-1999) and t₂ (2000-2018). (A) cloglog transformed output from full MaxEnt models for two time periods. Lighter yellow areas denote higher probabilities of occurrence. (B) Threshold maps of presence absence for the two time periods. Light areas represent predicted occurrence and dark gray represents predicted absence. Purple crosses represent actual occurrence data for each time period.

FIGURE 3

Figure 3. MaxEnt model predictions for predominant giant swallowtail butterfly host plants (Z. americanum, Z. clava-herculis and P. trifoliata) for T1 (1959-1999) and T2 (2000-2018). (A) cloglog transformed output from full MaxEnt models for each host plant across two time periods. Lighter yellow areas denote higher probabilities of occurrence. (B) Threshold maps of presence absence for the two time periods. Different colors (red, blue, and yellow) represent areas of predicted occurrence for each host plant and white represents predicted absence. Mixed colors indicate areas of overlap.

FIGURE 4

Figure 4. Kernel density estimates for modeled predicted presence of giant swallowtail butterfly P. cresphontes and host plants (Z. americanum, Z. clava-herculis and P. trifoliata) for T1 (1959-1999) and T2 (2000-2018). Red plots are from T1 (1959-1999) and blue plots are from T2 (2000-2018).

FIGURE 5

Figure 5. Predicted maximum northern occurence for P. cresphontes and host plants for T1 (1959-1999) and T2 (2000-2018). Dashed vertical lines represent the median value for each group. (A) P. cresphontes northern range limit between the two time periods. (B) Z. americanum northern range limit between the two time periods. (C) Northern range limit comparison for T1 (1959-1999) for P. cresphontes and Z. americanum. (D) Northern range limit comparison for T2 (2000-2018) for P. cresphontes and Z. americanum.

Papilio cresphontes Has Expanded Northward Due to Recent Climate Warming

Predictive maps generated from MaxEnt models clearly show a change in the distribution of P. cresphontes between T1 and T2, with a northward expansion since 2000 (Figures 2A,B). Kernel density estimate plots generated from threshold occurence predictions mirror this result (Figure 4), and highlight that different parts of P. cresphontes’ range match host plant use. Z. americanum closely matches P. cresphontes in the north, while the middle and southern part of the range is defined by the presence of Z. clava-herculis and P. trifoliata.

Host Plant Range Shifts

Overall, host plants (Z. americanum, Z. clava-herculis and P. trifoliata) demonstrated more complex changes in distribution between T1 and T2 compared to P. cresphontes (Figures 3A,B). Historically, the species were split latitudinally (with significant overlap) with Z. americanum occupying the northern part of the study area, P. trifoliata the middle, and Z. clava-herculis in the far south (Figure 3B). However, this pattern changes subtly in T2, with a range expansion of Z. americanum northward, but also westward to the boundary of our study area. Distribution changes in other host plants were more complex, with complicated range changes for P. trifoliata in the middle latitudes of the study area, and small range contraction of Z. clava-herculis to the south.

Northern Range Limits for P. cresphontes Have Shifted Northward and Closely Match Z. americanum

The northern range limit of P. cresphontes was significantly higher in T2 compared to T1 (t = −38.181, df = 560, p < 0.001; Figure 5A) where the median northern-most occurence for T2 (median = 46.1875 ± 0.675°) was 2.917° (∼324 km) higher in latitude than T1 (median = 43.2708 ± 1.692°). Z. americanum also demonstrated a significant (but small) northern range shift between T1 and T2 (t = −6.5717, df = 5510, p < 0.001; Figure 5B) where the median northern-most occurence for T2 (median = 45.5208 ± 0.914°) was 0.458° (∼51 km) higher in latitude than T1 (median = 45.0625. ± 1.667°). We also tested whether the northern range limits of P. cresphontes and Z. americanum differed from each other during each time period. In each time period, there was a significant difference between the northern range limits of P. cresphontes and Z. americanum (T1: t = −17.485, df = 550, p < 0.001; T2: t = 16.771, df = 551, p < 0.001). The difference between median butterfly and host plant northern range limits shrank from 1.75° (∼194 km) in T1 (with Z. americanum having a higher northern range limit) to 0.77° (∼85.47 km) in T2 (with P. cresphontes having a slightly higher median northern range limit; Figures 5B,C).

Climatic Variation in the Study Area Between T1 and T2

Overall, T2 had a higher mean annual temperature (9.45 ± 6.20° C) than T1 (8.67 ± 6.27°C) (t = −45.274, df = 534850, p < 0.001). Bioclim variables 10 and 11 [mean temperature of warmest quarter (breeding season) and mean temperature of the coldest quarter (pupal overwintering season)] had the biggest impacts on predicting P. cresphontes distribution, while variables 9 (mean temperature of driest quarter), 10 (mean temperature of warmest quarter) and 3 (isothermality) had the biggest impacts across both time periods for Z. americanum. Other host plants had multiple bioclim variables across time periods that impact distribution models (Figure 6). Variables that commonly had high permutation importance scores showed significant differences between T1 and T2 on average across our study area, with an overall trend of warmer patterns from 2000 to 2015 (T2) compared to 1959-1999 (T1) (Table 2).

FIGURE 6

Figure 6. Percent contribution of each of the 19 bioclim variables to final models for each species and time period. Percentages are computed from MaxEnt model training – as predictive gains increase, environmental factors contributing to feature generation are calculated and summarized in the final model. Common major contributors across many models include BIO9 (mean temperature of the driest quarter), BIO10 (mean temperature of the warmest quarter) and BIO11 (mean temperature of the coldest quarter).

TABLE 2

Table 2. Bioclimatic shifts in Bioclim variables between T1 (1959-1999) and T2 (2000-2015) that impact butterfly and host plant distributions.

Discussion

The determinants of species distributions have long been debated not just because they are essential in ecology and evolutionary biology, but also because where organisms are and where they will be is central to successful conservation and restoration practices in light of rapid climate change (Buckley et al., 2013; Gallagher et al., 2013; Robillard et al., 2015). Our study details a recent and rapid northward range expansion by P. cresphontes between 2000 and 2018 (Figure 1). We also model the distributions of the butterfly’s naturally occurring larval host plants, which, when combined with analysis of P. cresphontes range, result in different conclusions for the future distribution of this butterfly than if we had relied on abiotic variables alone (Figures 2, 3). Recent climatic shifts, particularly warmer, wetter temperatures during breeding season and warmer temperatures during pupal overwintering season, have allowed P. cresphontes to rapidly expand northward to now match or even surpass the slower moving northward range expansion of the northernmost host plant, Z. americanum, with further northward expansion of P. cresphontes now limited by host plant range, not climate (Figure 4). Our results highlight the importance of including biotic interactions (and interactions between herbivorous insects and host plants in particular) in examinations of range shifts and their speed, an idea often highlighted, (Urban et al., 2016) but infrequently implemented (Lemoine, 2015; Dilts et al., 2019; Svancara et al., 2019).

Poleward range shifts in herbivorous insects, particularly butterflies, have been documented for a number of species (Parmesan et al., 1999; Warren et al., 2001; Pöyry et al., 2009; Breed et al., 2012). Additionally, northward expansions of other butterfly species have been shown to have dramatic impacts on community composition through linked biotic interactions (Audusseau et al., 2017), which could be happening in this system as well, but would require further examination to determine. While studies demonstrating range shifts in multiple taxa provide valuable insights into the magnitude and direction of shifts for different taxa, gaps in knowledge remain (Pöyry et al., 2009). Namely, (1) how has warming acceleration affected recent range shifts during the last 10-15 years in poleward latitudes, and (2) how do abiotic and biotic factors interact to shape range shifts? Our study addresses both of these questions and provides a scalable, data acquisition and analytic pipeline by focusing on a single herbivore and multiple host plant species. We show an usually rapid northward range shift in this insect herbivore, P. cresphontes, over the last 18 years (predicted most northward occurrences differ by 2.917° of latitude (∼ 324 km) between T1 and T2, or a northward expansion of 180 km/decade) that is more than 27 times faster than the average of northward movement of global meta-analyses for plants, lichens, birds, mammals, insects, reptiles and amphibians, fish and marine organisms (Parmesan and Yohe, 2003) and over nine times faster than all butterfly species in Britain (Hickling et al., 2006). These observations are associated with warmer, wetter climate conditions during active flight times and overwintering. Our findings largely follow (Pöyry et al., 2009), who postulate that mobile species utilizing woody host plants like P. crespontes should exhibit large and fast range shifts northward, and that habitat availability and dispersal capacity largely determine success. We have laid the groundwork for one way to gather large amounts of data and analyze it as scale for future work across all butterfly and host plant species.

Interestingly, the northward incursion of P. cresphontes in northeastern North America is not a new phenomenon. Accounts detail movement into the region 145 years ago that lasted several decades (Scudder, 1889). In 1875, P. cresphontes were found in southern New England and by 1882 there are documented records just south of Montreal, Quebec. By the 1930s, the species had apparently retracted southward and were considered “extremely rare” in Massachusetts (Farquhar, 1934) and did not push northward into the region again until the last 8 years. Multiple long-term climate reconstructions (paired with historic instrument data) for the 145-year incursion period indicate a strong warming trend compared to the previous century (Marlon et al., 2016). However, this warming trend continues through the 1930s, so it is unclear which factors may have resulted in a retraction, though hydroclimatic reconstructions indicate an increase in drought in the northeastern United States over this time period, which likely had strong impacts on host-plant/nectar-plant distributions and quality through the range of P. cresphontes (Marlon et al., 2016), not to mention direct impacts on insect survival.

Our work also highlights the importance of including biotic interactions when predicting and projecting range shifts. Papillio cresphontes’ current northern range now closely matches the northernmost host plant (Z. americanum) (Figures 3, 5D) and this butterfly species is now limited by the ability of Z. americanum to expand its range northward. Because of the differences in life-history strategies, dispersal capabilities, reproductive outputs and environmental tolerances between insect and host plant, the northern expansion of P. cresphontes appears to now be largely curbed as the host plant is much more sessile and has much longer generation times. Though sightings of the winged adult stage of P. cresphontes will likely continue to be seen further north than the naturally occurring host plant range (Figure 5D), without a suitable host plant, further northward expansion seems unlikely but may be facilitated by recently documented P. cresphontes occurrences in horticultural settings. Papilio cresphontes lay eggs and larvae feed successfully on two non-native garden plants, garden rue (Ruta graveolens) and gas plant (Dictamnus albus). Common hoptree (P. trifoliata), is increasingly planted as an ornamental in the Northeast yet is a native species from central and southeastern North America. Although these exotics are not distributed uniformly across the region, dispersing P. cresphontes have an uncanny ability to find host plants in complex environments, perhaps further enabling them to expand their range in urban and suburban areas as abiotic conditions allow.

Data from community science sources continue to grow as platforms become more popular, and can provide tremendous boons to researchers across disciplines (Bonney et al., 2009, 2014; Dickinson et al., 2010), including those interested in creating species distribution models (Kéry et al., 2010; Yu et al., 2010). There has been debate about the quality and veracity of community science data, but recent work has demonstrated that citizen science initiatives can reliably produce research quality data though it often has similar biases to professionally-gathered data (Kosmala et al., 2016). Here, we use community science data sources supplemented by data from museum collections to generate species distribution models using the well-established MaxEnt modeling framework (Phillips and Dudík, 2008; Elith et al., 2011; Phillips et al., 2017), and advocate for continued development and use of community science data and its pairing with museum collection data in developing species distribution models in ecology and conservation.

Though we focused mostly on the distributional changes of P. cresphontes, there were also surprisingly large range shifts in host plant species (Figure 3). In contrast to the straightforward northward expansion of P. cresphontes, the distributional changes in host plants were more complex and nuanced. Z. americanum and P. trifoliata have both shifted northward between the two time periods in slightly different patterns (Figures 3, 4). While P. trifoliata appears to have shifted mostly northward (primarily gone from a large southern zone in T1), Z. americanum has undergone a northward and westward shift, and occupies areas that overlap with the range of P. trifoliata (Figure 3). The potential effects of this overlap on P. cresphontes (i.e., population dynamics, apparent competition, selection for oviposition behavior) are to our knowledge currently unknown, and warrants further examination in light of P. cresphontes westward expansion and previous work demonstrating significant within-population variation in oviposition behavior in Papilio (Thompson, 1988). Interestingly, mean temperature and annual temperature range (Bioclim variables 1 and 7) had the strongest impact in predicting the distribution of Z. americanum in T2, highlighting the impact that temperature may have in shaping and limiting current distribution. In contrast, the range of Z. clava-herculis appears to have contracted slightly in the southern United States. Compared to pre-2,000 distributions, available host plants to P. cresphontes are more widely distributed with greater overlap, but with notable gaps throughout portions of the southern United States. These complex distributional changes are likely driving part of the overall range shift northward for P. cresphontes (Figure 1B) and could also be potential drivers of speciation, and the evolution of specialization or host plant switching now and in the future (Descombes et al., 2016).

Conclusion

Multiple biotic interactions have evolved between insects and other species to create a wide variety of ecosystem services including herbivory and pollination (Losey and Vaughan, 2006). Anthropogenic climate change and habitat loss are creating a growing urgency for quantifying range size, understanding range boundaries, and assessing range shifts across insect species in order to preserve the integrity of future ecosystem function. Our work outlines the power of using increasingly abundant citizen science data, as well as the importance of including biotic interactions alongside environmental factors when developing analytical pipelines for biodiversity benchmarking studies. Future work should also incorporate climate change estimates into modeling efforts to project future distributions for both herbivores and host plants across many more butterfly and plant species. Incorporating both abiotic and biotic interactions in biodiversity benchmarking will provide a deeper, more nuanced understanding of temporal and spatial overlap among species, guiding conservation and management practices in a rapidly changing climate.

Statement of Data Archiving

Data and R scripts for all analyses are archived on Zenodo (https://doi.org/10.5281/zenodo.4476735).

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Zenodo doi: 10.5281/zenodo.4476735 (also in mansucript under Statement of Data Archiving).

Author Contributions

ML, KM, KP, DB, and JK conducted the project conception. JW carried out analyses (data acquisition, model building, statistical analyses, and visualization) with initial explorations of data and analyses from NC and support from RH. JW led manuscript preparation, with initial pieces in place from NC and DB. All authors supported in editing, commenting and adding material to the manuscript.

Funding

RH was supported in part by the National Science Foundation under grant no. 1910118.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank our reviewers and editors and Jeff Oliver for helpful comments and suggestions on the manuscript and analyses. Specifically, we would like to thank the Global Biodiversity Information Facility, the Maine Butterfly Atlas, the Maritime Canada Butterfly Atlas, the Massachusetts Butterfly Club, Butterflies and Moths of North America and eButterfly for the use of their data. We would also like to thank Scott Laurie and Ken-ichi Ueda at iNaturalist for their web platform, and all those relentless volunteers in various citizen science programs, past and present, who contributed these data for the greater good. This manuscript has been released as a pre-print at bioRxiv, (Wilson et al., 2019).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2021.579230/full#supplementary-material

Supplementary Figure 1 | Hyperparameter tuning and model evaluation of P. cresphontes distribution during the T1 period.

Supplementary Figure 2 | Hyperparameter tuning and model evaluation of P. cresphontes distribution during the T2 period.

Supplementary Figure 3 | Hyperparameter tuning and model evaluation of Z. americanum distribution during the T1 period.

Supplementary Figure 4 | Hyperparameter tuning and model evaluation of Z. americanum distribution during the T2 period.

Supplementary Figure 5 | Hyperparameter tuning and model evaluation of Z. clava-herculis distribution during the T1 period.

Supplementary Figure 6 | Hyperparameter tuning and model evaluation of Z. clava-herculis distribution during the T2 period.

Supplementary Figure 7 | Hyperparameter tuning and model evaluation of P. trifoliata distribution during the T1 period.

Supplementary Figure 8 | Hyperparameter tuning and model evaluation of P. trifoliata distribution during the T2 period.

Supplementary Table 1 | Bioclim variable definitions.

Footnotes

References

Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A., and Hegewisch, K. C. (2018). TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958-2015. Sci. Data 5:170191. doi: 10.1038/sdata.2017.191

PubMed Abstract | CrossRef Full Text | Google Scholar

Altieri, M. A., Nicholls, C. I., Henao, A., and Lana, M. A. (2015). Agroecology and the design of climate change-resilient farming systems. Agron. Sustain. Dev. 35, 869–890. doi: 10.1007/s13593-015-0285-2