Herbaria Reveal Herbivory and Pathogen Increases and Shifts in Senescence for Northeastern United States Maples Over 150 Years

Recent studies suggest climate-related delays in the timing of leaf coloration and abscission in maple trees but lack baseline data prior to the late 20th century. To better understand how autumn foliar phenology and late-season damage risks have changed for this genus over the past century, we evaluated 2,972 digitized herbaria specimens of red and sugar maple collected between 1826 and 2016 for the presence of leaves, autumn leaf coloration, and pathogen or herbivore damage. We found that the onset (first appearance) of colored leaves has shifted 0.26 days later each year, leading to a delay of more than a month in autumn phenology since 1880. We find that these shifts are related to precipitation regimes in both the fall and summer seasons and that more severe droughts are associated with higher probabilities of colored leaves. Moreover, we found that the probability of both herbivory and pathogen damage has increased significantly over the study period. In particular, we find a strong association between increasing summer drought conditions and increased probability of herbivory. Furthermore, the presence of foliar damage increased the probability of leaf coloration on herbaria specimens. However, the end-of-season abscission date (last appearance of leaves) was strongly associated with herbivory and climate in a contrary direction: Increasing yearly drought, higher fall temperatures, and the presence of herbivory were associated with earlier abscission. In fact, the last leaf dates for specimens with herbivory were nearly 2 weeks earlier than specimens without herbivore damage. Our study documents significant changes in maple senescence over the last 150 years and suggests that incorporating herbivory into models may improve our ability to predict forest responses to climate shifts.


INTRODUCTION
Consistently, scientists find that climate change has lengthened the growing season of deciduous forests by advancing spring leaf out (Badeck et al., 2004;Wolfe et al., 2005;Menzel et al., 2006;Polgar and Primack, 2011;Wolkovich et al., 2012;Ellwood et al., 2013). However, findings for autumn phenology have varied, with some suggesting temperature-related advancements in leaf senescence and abscission (Zani et al., 2020), while others report delays (Norby et al., 2003;Chen et al., 2020). The length of individual studies may contribute to different findings. Short term experiments show a delay of autumnal leaf senescence and abscission for some species, such as maple, in response to soil warming and increasing ambient temperatures (Norby et al., 2003;Wheeler et al., 2016). Remote sensing studies also document a later end to the vegetation growing season of approximately 5 days since the 1980s (Liu et al., 2016). However, remote sensing-derived vegetation data are not available prior to the late 20th century (National Research Council, 2008;Ustin and Gamon, 2010), limiting the utility of these methods to elucidate long-term trends. In addition, long-term ecological studies of phenology are still relatively rare, so current field studies may not be able to detect slowly developing environmental changes (Magnuson, 1990;Hobbie et al., 2003). Because of the limited temporal scope of studies, many questions remain about the responses of forest ecosystems, including how long climate stress has been influencing forest phenology and if the rate is increasing in current years.
Because a significant number of herbaria and museum specimens were collected before the intensification of anthropogenic climate change, such specimens can provide a baseline for investigating the impacts of environmental change on trends of plant phenology (Primack et al., 2004;Everill et al., 2014;Yost et al., 2018) and species interactions (Coleman and Seybold, 2011;Beauvais et al., 2017;. Preserved specimens provide an opportunity to look at a longer, more detailed time series of plant phenology than remotesensing data. For example, they may provide species-specific information missing from satellite measures that estimate the phenology of the canopy tree community as a whole. Studies using preserved specimens to investigate phenological changes are increasing as digitization enables access to large volumes of preserved specimens (Calinger et al., 2013;Davis et al., 2015;Besnard et al., 2018;Garretson, 2021). However, published studies focus on spring and summer flowering times and only a handful address autumn phenology (Forkner, 2013;Gallinat et al., 2015;Zani et al., 2020) and species interactions (Lavoie, 2013;Jones and Daehler, 2018;Meineke et al., 2018;Heberling and Isaac, 2019;Lang et al., 2019). Highly detailed studies that have examined leaf senescence have been limited to a few recent years (Panchen et al., 2015).
Herbaria specimens may confirm patterns of change but can also provide additional information missing from historical data by documenting insect and pathogen damage or plant morphological change. Thus, herbaria specimens can help in an overlooked area of phenology: the role of biotic factors and species interactions in augmenting or offsetting abiotic drivers. Insect damage, for instance, can cause greenfall or advance senescence (Forkner, 2013). Warming is generally expected to increase insect abundance, growth, and survival, particularly in temperate and boreal climates (Bale et al., 2002;Deutsch et al., 2008;Kingsolver et al., 2013) and, thus, also insect damage. However, at lower latitudes and in areas with higher average temperatures, summer temperatures may exceed insect thermal optima and reduce fitness or abundance (Deutsch et al., 2008). Furthermore, increasing summer and winter temperatures may have counteracting impacts on insect populations depending on the timing of insect life-history traits such as diapause, pupation, overwintering stage, or egg hatch. Winter temperatures have increased more than other seasons, which may be critical in changing insect survival and fitness (Bale et al., 2002;Elad and Pertot, 2014). Climate change studies of insects of forest trees have focused primarily on invasive species, outbreaking insects that can cause full-scale tree species loss, and on insect population dynamics. These studies rarely investigate insects as drivers of plant phenological processes. Like drought and high temperature, herbivory stress may alter plant physiology in ways that may impact autumn phenology. We do not have precise predictions because there are few experimental studies of this phenomenon (Belovsky and Slade, 2000;Forkner, 2013), but there is evidence from herbaria specimens that insect damage is increasing on plants .
Climate change is also likely to alter plant viral, fungal, and bacterial epidemics in a variety of ways. Studies show more intense pathogen outbreaks depending on pathogenspecific and host-specific characteristics (Ayres and Lombardero, 2000;Dale et al., 2001;Dukes et al., 2009;Jones and Barbetti, 2012;Jactel et al., 2012). Because host plants, insects, and pathogens will each be impacted differently by climate change, ecologists expect significant regional and species-specific changes in the type, intensity, and spatial distribution of plant damage (Shaw and Osborne, 2011;West et al., 2012;Chakraborty, 2013;Juroszek and von Tiedemann, 2013;Elad and Pertot, 2014;Garrett et al., 2016). Regardless of the type of damage, however, deciduous trees may prematurely abscise damaged leaves, change the timing of leaf maturation or bud or shoot production, or increase autumn pigmentation in response to damage (Williams-Linera and Baltazar, 2001;Forkner, 2013;Zvereva and Kozlov, 2014).
The impacts of interactions between abiotic and biotic forces on autumn phenology are likely to be complex. Difficulties in predicting separate impacts and their consequences complicate efforts to forecast changes in forest phenology. Specific trends in drought persistence and frequency are unclear (Dai et al., 2004;Huntington, 2006;Sheffield et al., 2012;Dai, 2013;Trenberth et al., 2014). Continued climate change will likely lead to increased drought (Cook et al., 2015;Touma et al., 2015) and, thus, to water stress in forest trees. Autumn coloration in deciduous trees may be delayed in response to summer drought stress (Xie et al., 2018b), while autumn drought leads to more species-specific responses. In one study, drought delayed peak coloration for A. rubrum but advanced autumn color in A. saccharum (Xie et al., 2018a). Drought and general changes in precipitation frequency and amount also influence mortality and morbidity of trees due to subsequent changes in herbivory and pathogen damage (Dale et al., 2001;Jactel et al., 2012;Marquis et al., 2019). Moderate water stress may benefit some herbivores, including leaf chewers or specialist leaf miners (Gely et al., 2020). Outcomes vary for particular pathogens, with some pathogens decreasing with drought and others increasing (Desprez-Loustau et al., 2006). Estimates of the impacts of drought on leaf senescence, the impacts of herbivores and pathogens on senescence, and the interactions of these factors are necessary to model future changes in tree autumn phenology.
To establish long-term patterns for both abiotic and biotic stressors on autumn phenology and determine whether insect and pathogen damage have changed with climate warming, we examined herbaria specimens of red and sugar maples collected over the past two centuries. Here we investigate associations of temperature, drought, herbivory, and pathogen damage on the timing of leaf abscission and senescence for maples in the eastern United States between 1826 and 2016. We documented changes in the leaf area removed and other types of damage on herbaria specimens and related these to abiotic factors over the collection dates to determine if herbivory and damage were associated with alterations in the magnitude and direction of shifts in autumn phenology.

Study Species
Acer rubrum and Acer saccharum (Sapindaceae) are ideal species to evaluate long term climate changes, herbivory, and pathogens using herbaria specimens because they (1) are widespread, abundant, and native to much of the eastern and central United States, (2) have unmistakable leaf shapes that reduce errors of identification and facilitate identification of herbivore damage, and (3) are well represented in herbaria collections in most months. Additionally, the tendency of these species to have red (anthocyanin and carotenoid) autumn coloration meant a higher likelihood of detecting this phenophase compared to other autumn colors, as dried specimens tend to brown, not produce red or orange color, after collection. Both species are medium to tall (20 to 40 meters) deciduous hardwoods with lifespans of 200 to 400 years (Loehle, 1988;Lorimer et al., 2004). Both have distinct, early season flowering and fruiting phenologies, producing flowers and samaras from March to April. Acer rubrum is one of the earliest flowering hardwoods and is more tolerant than A. saccharum of wet soils and high temperatures. The species show slight differences in autumn coloration timing, with A. saccharum senescing slightly later than A. rubrum (Archetti et al., 2013). Both show yellow to orange-red or bright red autumn color, with research suggesting phenotypic variation in the red intensity related to temperature (Xie et al., 2018a) and foliar nitrogen (Schaberg et al., 2003) for A. saccharum. Pathogens of A. rubrum and A. saccharum include tar spot fungi (Rhytisma spp.) and several fungi that cause anthracnose or other fungal rot (e.g., Colletotrichum, Discula, Kabatiella, and Phyllosticta). While the diversity of maples herbivores is lower than that of other dominant deciduous trees of the region, such as white oak (Quercus alba) (Marquis et al., 2019), several Lepidoptera of varying phenology cause notable chewing damage, including early season Tortricidae (Archips, Choristoneura, and Sparganothis) and Lasiocampidae (Malacosoma disstria Hübner), early and mid-season Geometridae and Noctuidae, and lateseason Noctuidae (Acronicta and Morrisonia), Notodontidae (Heterocampa and Symmerista), and Saturniidae (Dryocampa rubicunda Fabricius). However, of note for this study is that damage by one Hymenoptera, the maple eyespot gall [Acericecis ocellaris (Osten Sacken), Cecidomyiidae] was difficult to distinguish from fungal spot pathogens in older herbaria specimens.

Herbaria Data
We queried the Southeast Regional Network of Expertise and Collections (SERNEC; Gries et al., 2014;Zomlefer et al., 2017) for digitized specimens of A. rubrum and A. saccharum from the Central and North Eastern United States. Our query yielded 3,214 records distributed across 644 counties in Connecticut, Delaware, the District of Columbia, Illinois, Indiana, Kentucky, Maine, Maryland, Massachusetts, New Hampshire, New Jersey, New York, North Carolina, Ohio, Pennsylvania, Rhode Island, Tennessee, Vermont, Virginia, and West Virginia (Figure 1 Specimen dates ranged from 1826 to 2016. We then visually examined each record to verify the databased metadata using the specimen label data, focusing specifically on verifying the species, county of collection, and collection date. We excluded specimens without a month or year of collection, as well as specimens without a documented collection state or county of collection, from further analysis. We then georeferenced all points to the county centroid because only 17 specimens (0.5%) had exact GPS points at the time of data export. The number of records per county ranged from 1 to 41 specimens, with an average of 4.5 records per county (Figure 1). Additionally, we excluded specimens without visible pressed plants (e.g., those that included only seed packets or envelope preparations) or specimens with too poor an image quality to be characterized. Finally, we excluded specimens that mentioned being grown in a greenhouse or cultivated specimens to focus on wild specimens.
We visually annotated the remaining images for the presence or absence of leaves and colored leaves following the Nature's Notebook protocols for species-specific phenophase definitions and protocols described by Denny et al. (2014; Figure 1). These were developed for observations made on living plants, so phenophases are often excluded from analyses if the target plant part is dried or wilted. However, drying and wilting occur in the preservation of herbaria specimens, so we did not exclude these specimens in our scoring protocol. In our protocols, we annotated all evidence of red, yellow, or orange coloration as "colored leaves, " and any evidence of leaves as "leaves, " regardless of the collection date as the retention of leaves over winter could not be differentiated, and all efforts were taken not to bias phenophase assignment by the collection date. Described in Denny et al. (2014), the definition for colored leaves was: "One or more leaves show some of their typical late-season color (red, yellow, or orange) due to drought or other stresses. Do not include small spots of color due to minor leaf damage." The definition for leaves was: "One or more live, unfolded leaves are visible on the plant. A leaf is considered unfolded once its entire length has emerged from a breaking bud, stem node or growing stem tip, so that the leaf stalk (petiole) or leaf base is visible at its point of attachment to the stem." We also evaluated each specimen for the presence of leaf area removal caused by chewing damage by herbivores with mandibles, likely including FIGURE 3 | Scatter plots of the first day of year between the approximate summer and approximate fall solstice when a specimen was recorded with colored leaves (left) and the last day of year from the spring equinox to the spring equinox of the following year when a specimen was recorded with leaves (right). Points are colored by the species and sized by the number of specimens collected during that period. Lines are best fit linear regressions with a shaded area indicating the 95% confidence interval and are fit independently to each species.
primarily larval sawflies, butterflies, and moths (Hymenoptera and Lepidoptera), grasshoppers and allies (Orthoptera and Phasmid), and beetles (Coleoptera) (Figure 1). Leaf area removal was assessed because it is the most common type of herbivory easily recorded from herbaria specimens  and because it is the subject of most ecological studies on herbivory (Turcotte et al., 2014). Additionally, we annotated specimens for evidence of other pathogen damage, including bacterial, fungal, viral, or herbivore damage other than leaf area removal, but excluding damage due to pressing or tearing of leaves that may occur in the preservation of specimens.
We also excluded specimens from the same location, date, and collector with the same phenological and damage status to prevent overrepresentation bias based on collecting trips. Using a bias detection method similar to , we analyzed a subset of our data including only specimens from collectors that collected five or more specimens, a total of 62 collectors, and tested the effect of the collector on coloration, herbivory, and other damage using ANOVA. We found no evidence of collector bias in coloration collection patterns (F = 1.19, P = 0.153), but contrary to , we did find limited evidence of collector bias in the collection of damaged specimens. While changes to collection practices over time could confound assessments of foliar damage over time, we found this bias in collectors was not related to the mean year in which they collected or focused in earlier or later collectors, but rather was related to their season of collection. Specifically, collectors who differentially collected damaged specimens had later mean collection dates, and collectors who had differential collected undamaged specimens had earlier mean collection days of year. Regardless, we do not expect variation in collection preferences to vary by the variables of interest (e.g., we do not expect collectors in different localities or climates to show different biases), but as collectors are known to prefer collecting perfect specimens, our analyses likely show conservative or underestimates total foliar damage.

Climate Data
For specimens collected after 1895, we extracted mean annual and maximum fall (September to November) temperatures and precipitation at the county centroid from the Parameter-elevation Regressions on Independent Slopes Model (PRISM) Climate Data Explorer (PRISM Climate Group, 2020) for the year of collection. The annual temperatures included both daytime and nighttime data in mean temperature determinations. If the specimen was collected before September, the temperatures and precipitations were used from the fall of the prior year. Additionally, for each of the 644 counties included in the specimen dataset, 30-year normals (1981-2010) were extracted. We compared actual rainfall over the course of the year and specifically during the summer months to the yearly and summer rainfall normals by calculating The Standardized Precipitation Index (SPI; McKee et al., 1993). The SPI is the number of standard deviations by which the observed cumulative precipitation deviates from the climatological average. Positive values indicate wetter conditions compared to previous years, whereas negative values indicate drier years. As the baseline, we used countyspecific average monthly precipitation values over the 75 years from 1895-1970, as a comparison period of greater than 50 years is recommended (Guttman, 1999). As part of the PRISM export, we also extracted elevations for each of the locations. We Parameter estimates are reported ± standard error, along with variance inflation factors (VIF). The elevation is given in meters, temperatures in • C, and precipitation in millimeters. The SPI is the z-score of how much wetter or drier the time period was compared to the average precipitation during the summer at that location between . Negative values indicate significantly drier, while positive values indicate significantly wetter years. Reported coefficient Species: saccharum is for A. saccharum specimens relative to A. rubrum. *indicates <0.05, P > 0.01, **indicates <0.01, P > 0.001. *** indicates P < 0.001. excluded a total of 78 specimens collected before 1895 from climate-related analyses because PRISM does not have stable climate data available before 1895, and we excluded an additional 13 specimens that could not be georeferenced to county-level.

Statistical Analysis
We performed all analyses in R (version 4.0.2, R Core Team, 2016) with RStudio (version 1.2.5033, RStudio Team, 2020). Because phenology is a cyclical phenomenon that spans the end of the calendar year, evaluating phenophases at the end of the year linearly can lead to unintuitive results (Morellato et al., 2010;Pabon-Moreno et al., 2020). Because we are focused on autumn phenology, we shifted dates occurring before the approximate spring equinox (day 79, March 20th) to the end of the time series by adding 365, creating a new time series spanning from day of year 79 to 444. This modification allows for end dates of fall phenophases to extend into the following year without circular disruption (Figure 2).
To identify shifts in the timing of abscission and coloration, we analyzed the last day of year in each year in which any specimen showed the presence of leaves (a measure of end of the growing season) and the first day of year in each year with any specimen with presence of colored leaves (a measure of onset of senescence). We realize these first-and last-day phenophase metrics are finer-scale metrics than those typically used, such as peak or average phenophase. However, analyses suggest little difference in conclusions drawn from finer-scale stages in phenology models (Ellwood et al., 2019). Because immature leaves in both species can display red coloration, we omitted dates prior to the approximate summer solstice (day of year 171) and dates after the approximate winter solstice (day of year 355) to restrict onset dates to dates surrounding the autumn season.
Before testing more comprehensive models, we used simple linear regressions to identify shifts in the timing of the last date with recorded leaves over time and the first date with recorded colored leaves over time. Because presence/absence outcome variables are binary, we used logistic regressions to assess changes in the presence/absence of damage categories and colored leaves in specimen records over time. These less complex models provide baseline estimates of changes in these phenomena over time, prior to accounting for climate-related or spatial variability.
We used General Additive Models (GAMs) to evaluate the relationship between these dates and smoothed environmental conditions (the SPI of the collection year and the summer preceding collection; the fall precipitation and minimum and maximum temperatures), geospatial features (latitude, longitude, and elevation), the species, the year of collection, and damage (the presence of leaf area removal or other damage), and points were weighted by the number of specimens collected during that year. An additive model incorporates smooth functions of predictive variables, allowing us to reveal and estimate nonlinear relationships between predictors and response variables (Yee and Mitchell, 1991). This modeling approach also allows for the simultaneous modeling of latitude and longitude as a twodimensional smooth curve to account for spatial autocorrelation. In all models, temporal features (year in both models and day of year in the binary GAMs), elevation, and factors (species identity or the presence of damage) were modeled as parametric linear predictors. Location (as latitude and longitude) and climate variables were modeled as smoothed functions. We expect non-linear relationships over space, as there may be local adaptation, and with climate variables because extremely high and low levels of precipitation and temperature may be outside an organism's climatic niche. We tested for concurvity of the smoothed predictors with an estimated threshold of 0.8. These models were implemented using the mgcv R library (version 1.8-33; Wood, 2004Wood, , 2011Wood, , 2017Wood et al., 2016) and visualized using mgcViz (version 0.1.6; Fasiolo et al., 2020) and the tidyverse suite of packages (version 1.3.0; Wickham and RStudio, 2019).
We then explored the relationship between senescence and several variables -environmental conditions (the SPI of the collection year and the summer preceding collection; the fall precipitation and maximum temperatures), geospatial features (latitude, longitude, and elevation), the species, the year of specimen collection, the day of year of specimen collection, and damage (the presence of leaf area removal or other damage)using a binary GAM fit to the presence or absence of colored leaves on each herbaria specimen. We used a Locally Weighted Scatterplot Smoothing (LOWESS, Trexler and Travis, 1993) regression to visualize the probability of leaves on a specimen. However, neither a logistic regression nor binary GAM were computed because of a low number of specimens without leaves and an extended period during the year with a 0% chance of specimens without leaves (Figure 2). Similarly, we fit two binary GAMs to the presence or absence of (1) leaf area removal and (2) other damage on herbaria specimen to explore the relationship between damage and geospatial features (longitude and elevation), the species identity, the year of collection, the day of year of collection (as the probability of damage increases the longer a leaf is on a tree), environmental conditions (the SPI of the collection year and the summer preceding collection, the annual total precipitation, and annual mean temperature), and the presence of the non-target damage category. In contrast to the model for leaf senescence, because damage can occur before the fall season, we used the mean yearly temperature and the total annual precipitation instead of fall-specific climate. In these models, location (as latitude and longitude), elevation, and climate variables were modeled as smooth functions. For these classification models (the binary GAMs), we report the area under the receiver operating characteristic curve (AUC), which is a measure of classification performance that ranges from 0, indicating a model with 100% incorrect prediction, and 1, a perfect classification model, and we calculated this metric using the pROC R library (version 1.17.0.1; Robin et al., 2011).

RESULTS
A total of 2,972 herbaria records were included in our resulting dataset (Figure 1). The two study species were approximately equally represented with 55% (n = 1,640) A. rubrum and 45% (n = 1,337) A. saccharum. The majority of specimens (70%, n = 2,100, Figure 2) included preserved green leaves, while the preservation of colored leaves was less frequent (15%, n = 447, numbers do not sum to 100% because some specimens contained only stems with buds, Figure 2). As would be expected, the probability of leaves and colored leaves varied over the yearly cycle (Figure 2) and showed seasonal patterns of occurrence. Leaf area removal was documented on 52% (n = 1,568) of specimens, while 35% (n = 1,063) of specimens had evidence of other damage (fungal, viral, bacterial, or gall spots). The resulting data are available in the Environmental Data Initiative repository under identifier edi.713 (Garretson and Forkner, 2021).

Shifts in Phenophase Timing
The results of a simple linear regression of the year of collection with the timing of leaf-offset (last day of year with leaves) showed a significant change over time with an estimated delay of 0.43 ± 0.09 days per year (P < 0.0001), leading to a total estimated delay of approximately 84 days between 1826 and 2016 (Figure 3). However, the fit of the simple model was low, with R 2 = 0.09, justifying the exploration of relationships with climate and damage features. The General Additive Model with expanded predictors explained 43.4% of the deviance in the last observation date of leaves (Table 1 and Supplementary  Figure 1). We observed non-linear relationships with the location of collection (edf = 19.16, F = 3.44, P < 0.0001). In particular, we observed earlier last leaf dates in more northern areas but delays in the central Ohio region (Figure 4A). Higher elevation was associated with earlier dates of last leaves, with each meter associated with 0.14 ± 0.04 day earlier last leaf dates (t = −3.92 P = 0.0001, Figure 4B). Leaf area removal was associated with an earlier last observation of leaves date. Specimens with damage were on average 19.26 ± 8.11 days earlier than specimens without damage (t = −2.37, P = 0.02, Figure 4C). Additionally, we observe relationships between last leaf dates and the yearly SPI (edf = 1.00, F = 9.15, P = 0.003) and maximum fall temperatures (edf = 2.46, F = 3.10, P = 0.026). Non-linearity was particularly evident in the relationship with the maximum fall temperature: at temperatures higher than approximately 25 • C, there was no marked change from the average dates, but as the maximum temperature increased from 17.8 to 25 • C, there was a negative relationship between the temperature and the last leaves date, with the last leaves date moving closer to the average date ( Figure 4D). Increasing drought was associated with earlier last leaf dates, while increasing wetness relative to historical averages was associated with later last leaf dates ( Figure 4E). These differences were quite large, with an SPI of −2 (severe drought) associated with last leaf dates nearly 30 days earlier than average and an SPI of 3 (extremely wet) associated with last leaves dates nearly 40 days later than the average (Figure 4E).
A simple linear regression of the year of collection with coloration onset (the first day of year with colored leaves) also revealed a significant delay of approximately 0.21 ± 0.89 days per year, suggesting a total delay of approximately 26 days between 1890 and 2015 (Figures 2, 3). As with the model of abscission, the model fit with only the year of collection as an explanatory factor was very low, with an R 2 of only 0.03, again justifying more complex models. Comparably, the model fit of the GAM with expanded predictors for the first appearance of colored leaves explained 45.5% of the deviance (Table 1 and Supplementary Figure 2). There was evidence of non-linear relationships with the summer SPI (edf = 6.59, F = 2.57, P = 0.01, Figure 5A). In years where rainfall was higher than average, coloration was delayed, until rainfall exceeded two standard deviations from the mean and coloration occurred much earlier. In years when rainfall was lower than average, coloration onset was later ( Figure 5A). The first appearance of colored leaves in autumn was delayed by 0.26 ± 0.12 days per year over the study period (t = 2.22, P = 0.029, Figure 3). Importantly, including additional predictors suggested that the simple model provided an underestimate of the total shifts in coloration timing. Finally, as elevation increased, coloration was advanced by 0.06 ± 0.02 days (t = −2.83, P = 0.006, Figure 5B).

Probability of Coloration
A simple logistic regression showed that the likelihood of coloration on collected specimens did not significantly vary over time (P = 0.3). However, GAMs revealed that the probability that a specimen possessed colored leaves was related to the summer SPI (edf = 1.29, χ2 = 4.29, P = 0.048, Figure 6A), following the trends seen in the timing of coloration ( Table 2 and Supplementary Figure 3). Because coloration is a phenophase that primarily occurs in the autumn, as expected, there was a significant positive relationship between the adjusted ordinal day of year and the probability that a specimen contains colored leaves (β = 0.008 ± 0.001, z = 10.11, P < 0.0001, Figure 6B). Additionally, we find a strong relationship to damage metrics (Table 2 and Figures 6C,D). Both leaf area removal The elevation is given in meters, temperatures in • C, and precipitation in millimeters. The SPI is the z-score of how much wetter or drier the time period was compared to the average precipitation during the summer at that location between 1895-1970.
(β = 0.36 ± 0.13, z = 2.84, P = 0.005, Figure 6C) and other damage (β = 0.27 ± 0.12, z = 2.23, P = 0.026, Figure 6D) were associated with an increase in the probability that a specimen included colored leaves. Additionally, relative to A. rubrum, there was a lower probability that A. saccharum specimens contained colored leaves on a given date (β = −0.62 ± 0.12, z = −5.18, P < 0.0001, Figure 6E). The area under the curve (AUC) for this model was 0.72, meaning this model performed significantly better than chance in predicting the coloration status of any given specimen.

Probability of Damage
The results of a simple logistic regression using only year as a covariate show that without accounting for climate variation, the probability of both damage categories appearing on a specimen increased significantly over the study period (Leaf area removed: β = 0.004 ± 0.001, P < 0.0001; Other damage: β = 0.008 ± 0.001, P < 0.0001, Figure 7). With an expanded model including spatial and climate variation ( Table 2 and Supplementary Figures 4, 5), we find that increasing herbivory is strongly related to summer drought and annual temperatures, which drive the interannual variation over the past 100 years (Figure 8). Conversely, the probability of other pathogen damage retained the year as a predictive term, indicating the selected climate variables do not account for the interannual variation and observed increasing trends. This finding suggests factors not included in our model but increasing over time (e.g., urbanization, pollution, or introduced species proliferation) may play a prominent role in the drivers of pathogen damage increases (Figure 9). The predictive capabilities of these models were relatively high, with an AUC of 0.63 for the probability of herbivory and an AUC of 0.77 for the probability of other damage. Thus, the models performed significantly better than chance at predicting the damage status of the observations. As expected, damage was also positively associated with the day of year: as the year progressed, the probability of damage increased (Leaf area removed: β = 0.006 ± 0.001, z = 10.55, P < 0.0001, Figure 8B; Other damage: β = 0.007 ± 0.001, z = 10.33, P < 0.0001, Figure 9B). For both categories of damage, the presence of the alternative damage on the specimen significantly increased the probability of the other damage category, which suggests that herbivory increases susceptibility to other pathogen damage or vice versa (Leaf area removed: β = 1.33 ± 0.95, z = 13.99, P < 0.0001, Figure 8C; Other damage: β = 1.36 ± 0.95, z = 14.39, P < 0.0001, Figure 9C). The probability of leaf area removal varied spatially (edf = 15.88, χ2 = 94.34, P < 0.0001, Figure 8A). The pattern of spatial variation in leaf area removal showed lower probabilities of herbivory in more northern and coastal regions ( Figure 8A). The probability of herbivory was also associated with the summer SPI (edf = 1.0, χ2 = 4.09, P = 0.043, Figure 8D) and the annual mean temperatures (edf = 1.0, χ2 = 5.13, P = 0.024, Figure 8E), with increasing drought associated with a higher probability of herbivory and increasing wetness associated with decreased probability of herbivory. Lower average annual temperatures were associated with a higher probability of damage, but higher annual temperatures were associated with decreased probability of herbivory. Unlike herbivory, other damage did not display strong associations with climate variation. However, there were patterns of spatial variation (edf = 17.45, χ2 = 64.45, P < 0.0001), with northern and inland regions showing higher probabilities of other damage (Figure 9A).

DISCUSSION
Our analyses of red and sugar maple phenology suggest that autumn coloration has been delayed by approximately 0.26 days/year since 1826, leading to a delay of more than a month in recent decades (Figures 2, 3). We find that the probability of coloration and the timing of fall color is related to drought severity, with more severe droughts increasing the probability of colored leaves ( Figure 6A). Additionally, we find signals of increasing herbivory and pathogen damage recorded on herbarium specimens (Figure 7). These changes in foliar damage were associated with both the probability of coloration and the timing of maple leaf abscission. Foliar damage, both herbivory and pathogen damage, was associated with a significantly higher probability that a specimen showed autumn coloration ( Table 2). Herbivory was also associated with a nearly three-week earlier abscission date compared to undamaged specimens. Together, these results indicate that foliar damage by insect herbivores and pathogens is a critical biotic determinant of the timing and expression of autumnal phenophases, and that in some cases it may counteract climate-driven delays.
We find that drought and precipitation regimes play a vital role in the timing of abscission and autumn coloration in maples. We find that more severe drought was associated with significantly earlier last leaf observation dates. Specifically, an SPI two standard deviations drier than historical averages was associated with a last leaf date 25 days earlier than average ( Figure 4E). Wetter conditions significantly delayed the last observation date, with conditions three standard deviations wetter than historical averages associated with last leaf dates more than a month later ( Figure 4E). In conjunction with a positive but non-linear relationship with temperature, this could indicate an important role of evapotranspiration in modifying abscission dates and lengthening the growing season. Some previous studies found that daytime warming in autumn was associated with earlier abscission because higher daytime evapotranspiration without a corresponding increase in precipitation increases tree water stress (Chen et al., 2020). The first observation of colored leaves also showed a strong association with drought ( Figure 5A), as did the probability that a specimen contained colored leaves ( Figure 6A). Global climate change is projected to increase the variance of precipitation, increasing the risk of both floods and droughts (Meehl and Tebaldi, 2004;Trenberth et al., 2014;Hubbart et al., 2016). Under future climate scenarios, the combined effects of temperature, drought, wetness, and leaf damage may stress trees and negatively affect forest health and phenological timing.
Our finding of a 19-day advance in the last observation day of leaves when specimens were damaged by foliar herbivores is consistent with experimental manipulations of A. rubrum that showed earlier abscission on trees with high levels of simulated herbivory (Forkner, 2013) and with previous studies demonstrating shorter leaf lifespans in herbivore-damaged leaves (Risley and Crossley, 1988). The probability of leaf coloration was also increased with the probability of both herbivory and other damage ( Table 2). Red leaf coloration has been proposed as a signal to protect against excessive herbivory, particularly in immature leaves (Karageorgou and Manetas, 2006;Karageorgou et al., 2008;Green et al., 2015), and some red coloration may be related to fungal pathogens (Coley and Aide, 1989). Experimentally, red maples have been shown to prematurely color damaged leaves (Forkner, 2013). Damage may, therefore,  be working in opposition to climate-related delays, and increased herbivory over the past century could offset observed extensions of the growing season. Because abscission was earlier as drought and herbivory increased, these stressors may reduce the predicted lengthening of the growing season based on models using remotely sensed data.
Changes in herbivory probabilities were associated with drought conditions and temperature, with more severe drought associated with higher probabilities of herbivory damage and higher temperatures associated with lower probabilities of herbivory damage (Figure 8D). The direction and magnitude of temperature impacts on insect herbivores are known to vary by species (Lemoine et al., 2014). Increasing temperatures are generally expected to increase insect survival (Meineke et al., 2013;Dale and Frank, 2017), but extremely high temperatures can decrease insect survival and increase developmental time (Rocha et al., 2017). Both insect and pathogen damage increased over the time period of the study, with leaf area removal associated with an increase in the probability of other pathogen damage on a specimen and vice versa ( Table 2). Experimental studies have demonstrated that some insect herbivores can be attracted to diseased host plants if infection increases the nutritional quality of the plant or alters chemical defenses (Stout et al., 2006;Srisakrapikoop et al., 2020). Herbivory also causes wounding of plant foliar tissue, which can introduce microbial pathogens or cause further damage (Savatin et al., 2014).
Increases in the probability of pathogen damage over time were not solely explained by climate variability, suggesting that other changing forest features over the past two centuries may be responsible. This might include the introduction of invasive species, land-use changes, and urbanization (Meineke et al., 2013(Meineke et al., , 2014Dale and Frank, 2017). However, because of the observed relationship between herbivory and pathogen damage, pathogen damage may show similar patterns with drought and temperature in more detailed studies. Wetter years are associated with many fungal plant pathogens (Desprez-Loustau et al., 2007;Giauque and Hawkes, 2016;Hubbart et al., 2016), but other fungal, viral, and bacterial diseases can be more successful in attacking stressed plants during droughts (Dale et al., 2001;Desprez-Loustau et al., 2007;Jactel et al., 2012). Therefore, future studies may benefit from explicitly investigating the identities of pathogens causing damage, as fungal, viral, and bacterial diseases may have different climate relationships. Additionally, explaining the observed spatial patterns is challenging because they may be confounded by urbanization, deforestation, changes in deciduous tree species composition, and invasive species introductions over the time span of our analysis. Accordingly, we further suggest that future studies using herbaria specimens incorporate measures of forest management and urbanization over time in models of phenological change to improve our understanding of spatial variability in phenology and foliar damage over time.
In summary, our results imply that increasing damage may mask shifts in the timing of coloration onset and may partially offset end-of-season delays in the future. Although sample sizes were relatively small in this study, as digitization efforts increase the availability of specimens, our ability to detect relationships between climate and phenological timing using herbaria will improve (Beaman and Cellinese, 2012;Soltis, 2017;Giraud et al., 2018). Clearly, heat, water, insect, and pathogen stressors collectively impact autumn phenology, and models of future forest growth may need to incorporate all three for improved predictions. Our study indicates that under some conditions, later coloration and earlier abscission may mean shorter autumns overall.

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: https://portal. edirepository.org/nis/mapbrowse?scope=edi&identifier=713.

AUTHOR CONTRIBUTIONS
AG and RF contributed to the conceptualization, methodology of the article, and analyzed and annotated herbaria specimens. AG curated and analyzed the data and wrote the original draft. RF administered the project, provided supervision, and provided additional writing, review, and editing. Both authors contributed to the article and approved the submitted version. We thank Lorelei D. Crerar and Rebecca Dikow, as well as Rachel A. Keuler and Alexandra Duffy, for comments on the manuscript and analyses. We also thank Andrea Weeks, Director of the Ted R. Bradley Herbarium at GMU, for advice when starting the project. Additionally, we acknowledge Kadiera S. Ingram for her initial assistance. Finally, AG thanks Geraldine Grant and the Biology Research Semester participants for their comments on the development of the project.