Phylogenetic patterns of species loss in Thoreau's woods are driven by climate change The Harvard community has made this article openly available. Please share how this access benefits you. Your story matters Citation Willis, C. G., B. Ruhfel, R. B. Primack, A. J. Miller-Rushing, and C. C. Davis. 2008. Phylogenetic Patterns of Species Loss in Thoreau’s Woods Are Driven by Climate Change. Proceedings of the National Academy of Sciences 105, no. 44: 17029–17033. doi:10.1073/ pnas.0806446105. Published Version doi:10.1073/pnas.0806446105 Citable link http://nrs.harvard.edu/urn-3:HUL.InstRepos:29362004 Terms of Use This article was downloaded from Harvard University’s DASH repository, and is made available under the terms and conditions applicable to Other Posted Material, as set forth at http:// nrs.harvard.edu/urn-3:HUL.InstRepos:dash.current.terms-of- use#LAA http://osc.hul.harvard.edu/dash/open-access-feedback?handle=&title=Phylogenetic%20patterns%20of%20species%20loss%20in%20Thoreau's%20woods%20are%20driven%20by%20climate%20change&community=1/1&collection=1/2&owningCollection1/2&harvardAuthors=c30265921a2698d382919ae653f63381&departmentOrganismic%20and%20Evolutionary%20Biology http://nrs.harvard.edu/urn-3:HUL.InstRepos:29362004 http://nrs.harvard.edu/urn-3:HUL.InstRepos:dash.current.terms-of-use#LAA http://nrs.harvard.edu/urn-3:HUL.InstRepos:dash.current.terms-of-use#LAA http://nrs.harvard.edu/urn-3:HUL.InstRepos:dash.current.terms-of-use#LAA Phylogenetic patterns of species loss in Thoreau’s woods are driven by climate change Charles G. Willisa, Brad Ruhfela, Richard B. Primackb, Abraham J. Miller-Rushingb, and Charles C. Davisa,1 aDepartment of Organismic and Evolutionary Biology, Harvard University Herbaria, 22 Divinity Avenue, Cambridge, MA 02138; and bDepartment of Biology, Boston University, 5 Cummington Street, Boston, MA 02215 Edited by Michael J. Donoghue, Yale University, New Haven, CT, and approved September 19, 2008 (received for review July 3, 2008) Climate change has led to major changes in the phenology (the timing of seasonal activities, such as flowering) of some species but not others. The extent to which flowering-time response to tem- perature is shared among closely related species might have important consequences for community-wide patterns of species loss under rapid climate change. Henry David Thoreau initiated a dataset of the Concord, Massachusetts, flora that spans !150 years and provides information on changes in species abundance and flowering time. When these data are analyzed in a phylogenetic context, they indicate that change in abundance is strongly corre- lated with flowering-time response. Species that do not respond to temperature have decreased greatly in abundance, and include among others anemones and buttercups [Ranunculaceae pro parte (p.p.)], asters and campanulas (Asterales), bluets (Rubiaceae p.p.), bladderworts (Lentibulariaceae), dogwoods (Cornaceae), lilies (Lil- iales), mints (Lamiaceae p.p.), orchids (Orchidaceae), roses (Rosa- ceae p.p.), saxifrages (Saxifragales), and violets (Malpighiales). Because flowering-time response traits are shared among closely related species, our findings suggest that climate change has affected and will likely continue to shape the phylogenetically biased pattern of species loss in Thoreau’s woods. conservation ! extinction ! phenology ! phylogenetic conservatism ! phylogeny The impact of climate change on species and communities hasbeen well documented. Arctic forests are shifting poleward and alpine tree lines are shifting upward (1–3); spring f lowering time is advancing rapidly (4 –7); pest outbreaks are spreading (8); and numerous species are declining in abundance and risk extinction (9). However, despite these generalized trends, spe- cies vary dramatically in their responses to climate change. For example, although the spring f lowering times of many temperate plants are advancing, some are not changing and others are f lowering later in the season (5, 10, 11). Understanding the evolutionary (i.e., phylogenetic) history of traits that are inf lu- enced by climate (e.g., f lowering phenology) has been an un- derexplored area of climate change biology, despite the fact that it could prove especially useful in predicting how species and communities will respond to future climate change. Closely related species often share similar traits, a pattern known as phylogenetic conservatism (12–16, 17). If closely related species share similar traits that make them more susceptible to climate change (14, 17), species loss may not be random or uniform, but rather biased against certain lineages in the Tree of Life (i.e., phylogenetic selectivity; see ref. 18). However, a deeper inquiry into these patterns has been hampered largely because adequate datasets documenting community-wide responses to climate change are exceedingly rare. During the mid-19th century, the naturalist and conservationist Henry David Thoreau spent decades exploring the temperate fields, wetlands, and deciduous forests of Concord, Massachu- setts, in the northeastern United States. He wrote extensively about the natural history of the area (19) and kept meticulous notes on plant species occurrences and f lowering times (11, 20). Several botanists have since resurveyed the Concord area, thus providing a unique community-level perspective on changes in its f loristic composition and f lowering times during the past !150 years (11, 20). Despite the fact that !60% of all natural areas in Concord are undeveloped or have remained well protected, a striking number of species have become locally extinct: 27% of the species documented by Thoreau have been lost, and 36% exist in such low population abundances that their extirpation may be imminent (20). Also, the species that have been lost are overly represented in particular plant families (20), suggesting that extinction risk may be phylogenetically biased. Although habitat loss due to succession and development (e.g., loss of wetlands, abandonment of farms, reforestation, and construction of homes and roads) has contributed to decreases in abundance for some species in Thoreau’s Concord (20), climate change may also help to explain the seemingly nonran- dom pattern of species loss among certain plant groups. It has been shown recently (11) that the mean annual temperature in the Concord area has risen by 2.4 °C over the past !100 years and that this temperature change is associated with shifts in f lowering time: species are now f lowering an average of 7 days earlier than in Thoreau’s time. Along with changes in f lowering phenology, species range is likely to be inf luenced by climate change (21). Thus, the Concord surveys provide a unique opportunity to examine the extent to which changes in abun- dance may be correlated with these climatologically sensitive traits. Also, by incorporating phylogenetic history into our analyses, we can test whether species that share similar traits are closely related (i.e., phylogenetic conservatism), and to what extent these traits correlate with decreases in abundance. Such findings could identify groups of closely related species that are at higher risk of extinction (18, 22). The data for the 473 species we analyzed were collected by Thoreau (1852–1858), Hosmer (1878, 1888 –1902), and Miller- Rushing and Primack (2003–2007) (see Materials and Methods; see refs. 20 and 23). Scorings include information on changes in species abundance, species habitat, and 2 separate measures of f lowering-time response to temperature (i.e., the ability of species f lowering time to track short-term seasonal temperature changes, and the shift in species f lowering time over long-term intervals). We further scored the current mean latitudinal range and native/introduced status of each species. We constructed a composite phylogeny of all species to test for: (i) the phylogenetic conservatism of each trait, and (ii) correlations between these traits and change in abundance when accounting for phylogeny. Author contributions: C.G.W. and C.C.D. designed research; C.G.W., B.R., R.B.P., A.J.M.-R., and C.C.D. performed research; C.G.W., B.R., and C.C.D. analyzed data; and C.G.W., B.R., R.B.P., A.J.M.-R., and C.C.D. wrote the paper. The authors declare no conflict of interest. This article is a PNAS Direct Submission. Freely available online through the PNAS open access option. 1To whom correspondence should be addressed. E-mail: cdavis@oeb.harvard.edu. This article contains supporting information online at www.pnas.org/cgi/content/full/ 0806446105/DCSupplemental. © 2008 by The National Academy of Sciences of the USA www.pnas.org"cgi"doi"10.1073"pnas.0806446105 PNAS ! November 4, 2008 ! vol. 105 ! no. 44 ! 17029 –17033 EV O LU TI O N http://www.pnas.org/cgi/content/full/0806446105/DCSupplemental http://www.pnas.org/cgi/content/full/0806446105/DCSupplemental Results and Discussion Our results (Fig. 1 and Table 1) indicate that change in abundance and f lowering-time response traits were phyloge- netically conser ved, which indicates that species evolutionar y histor y is important to understanding community response to climate change. Species that are declining in abundance are more closely related than expected by chance. Similarly, species that exhibit similar f lowering-time responses to tem- perature are more closely related than expected by chance. In contrast, latitudinal range was not phylogenetically conser ved Fig. 1. Composite phylogeny of 429 flowering plant species from the Concord flora depicting changes in abundance from 1900 to 2007. Change in abundance ranged on an integer scale from "5 to #4, and was calculated as the difference in abundance for each taxon in 1900 and 2007 based on 7 abundance categories (0 to 6; see Materials and Methods). Branch color indicates parsimony character state reconstruction of change in abundance. For simplicity, we have indicated this reconstruction by using 4 colors: red (major decline, "5 to "3), pink (moderate decline, "2), gray (little to no change, "1 to #1), and blue (increase, #2 to #4). For the complete character reconstruction and taxon labels see Fig. S1. Average decline in abundance was calculated for all internal nodes as the mean change in abundance of descendant nodes weighted with branch length information ascertained from divergence time estimates. An average decline of 2.5 or greater corresponds to a decline in abundance of 50% or greater, based on our most conservative scoring using 6 abundance categories (0 to 5; see Materials in Methods). Clades exhibiting these major declines are indicated with black dots. Each of the most inclusive clades exhibiting these declines are indicated in pink and referenced numerically to their clade name. Subclades in major decline that are nested within more widely recognized clades are labeled with the more familiar name followed by pro parte (p.p.). These clades include some of the most charismatic wildflower species in New England, such as anemones and buttercups (Ranunculaceae p.p.), asters, campanulas, goldenrods, pussytoes, and thistles (Asterales), bedstraws and bluets (Rubiaceae p.p.), bladderworts (Lentibulariaceae), dogwoods (Cornaceae), lilies (Liliales), louseworts and Indian paintbrushes (Orobanchaceae), mints (Lamiaceae p.p.), orchids (Orchidaceae), primroses (Onograceae p.p.), roses (Rosaceae p.p.), saxifrages (Saxifragales), Indian pipes (Ericales p.p.), and St. John’s worts and violets (Malpighiales). Table 1. Statistical tests of phylogenetic conservatism and trait correlations with change in abundance Trait Phylogenetic conservatism Trait correlation Model 1 Model 2 Model 3 n Observed rank n Estimate n Estimate n Estimate Flowering time tracking of seasonal temperature 175 19 ** 175 "0.48 * 166 "0.62 * 140 "1.00 *** Shift in flowering time 1850–1900 319 2 *** 319 "0.02 *** 311 "0.01 * 140 0.03 *** Shift in flowering time 1900–2006 303 2,120 — 303 0.04 *** 296 0.03 *** 140 0.02 *** Shift in flowering time 1850–2006 271 340 † 271 0.04 *** 253 0.03 *** 140 — — Mean latitudinal range 414 3,705 414 "0.10 *** 362 "0.08 *** 140 "0.09 *** Change in abundance 1900–2006 429 1 *** — — — — — — — — — Tests used a phylogeny with branch lengths adjusted for time. The significance of phylogenetic conservatism was tested by comparing the rank of the observed standard deviation (SD) of descendent trait means to a null model based on 9,999 random iterations of trait distributions across the composite phylogeny. The observed rank is compared with a 2-tail test of significance, i.e., an observed rank of 250 equals a P value of 0.05. Trait correlations were tested by using the comparative methods of generalized estimating equations (GEE). Estimates describe the direction and magnitude of the correlation (e.g., a negative estimate $"0.1% of mean latitude with change in abundance suggests that species from more southerly latitudes are increasing in abundance). Model 1 (univariate model), correlation of change in abundance with each trait; Model 2 (multivariate model), correlation of change in abundance with each trait and habitat, abundance (ca. 1900), flowering season, and native/introduced status as covariates; Model 3 (multivariate model), correlation of change in abundance with all traits and habitat, abundance (ca. 1900), flowering season, and native/introduced status as covariates (shift in flowering-time response 1850 –2006 was excluded due to its high correlation with the other flowering-time shift traits). †, P & 0.1; *, P & 0.05; **, P & 0.01; ***, P & 0.001; n & sample size. 17030 ! www.pnas.org"cgi"doi"10.1073"pnas.0806446105 Willis et al. http://www.pnas.org/cgi/data/0806446105/DCSupplemental/Supplemental_PDF#nameddest=SF1 (i.e., phylogeny is not important in explaining the latitudinal distribution of species). The ability of species to track seasonal temperatures was correlated with changes in abundance: species whose f lowering time does not track seasonal temperature have greatly declined in abundance over the past !100 years. Similarly, shifts in species f lowering time across all 3 long-term time intervals (1850 –1900, 1900 –2006, and 1850 –2006) were correlated with change in abundance: species that are not f lowering earlier have declined in abundance. Last, species range was correlated with change in abundance: more northerly species have decreased in abundance in relation to southerly species. Our results are robust: (i) when controlling for multiple variables that may additionally affect decline in abundance [i.e., initial abundance, habitat, native/ introduced status, and f lowering season (date of first f lowering); see Table 1], (ii) to branch length information [supporting information (SI) Table S1], and (iii) to phylogenetic uncertainty (Table S2). These results demonstrate that there is a phylogenetically selective pattern of change in abundance. Decreases in abun- dance have been disproportionately high in certain clades, including asters, bladderworts, buttercups, dogwoods, lilies, louseworts, mints, orchids, saxifrages, and violets (see Fig. 1). This result confirms previous f loristic studies across similar time spans demonstrating that the risk of plant extinction (i.e., occurring in low abundance; see ref. 24) is taxonomically (20, 25–27) and phylogenetically (28) shared among close relatives. However, to our knowledge our study is the first to report that the phylogenetic selectivity of extinction risk is correlated with traits directly inf luenced by climate change. Species whose f lowering times are not responsive to changes in temperature are decreasing in abundance. Most strikingly, species with the ability to track short-term seasonal temperature variation have fared significantly better under recent warming trends. In addition, species whose f lowering times have shifted to be earlier in the year over the long-term have also fared significantly better under recent warming trends. Based on our regression estimates (Table 1), change in abundance over the last !100 years is greatest when assessed against the ability of species to track short-term sea- sonal temperature versus long-term f lowering shifts. Thus, the association between f lowering-time tracking and change in abundance is a better estimator of species response to rising temperatures. Interestingly, these 2 f lowering-time response traits are significantly, but weakly correlated. This weak corre- lation raises the possibility of different mechanisms of pheno- logical response to climate change (e.g., plasticity, adaptation; see refs. 29 and 30). Alternatively, confounding factors such as changes in population size may affect estimates of long-term shifts in first f lowering dates, but would be less likely to inf luence estimates of tracking climate change over the short term (31). Asynchronous phenological responses resulting from rapid climate change can have negative fitness effects on organisms, leading to dramatic declines in population sizes or local extinc- tion (32). Selection on f lowering phenology may be direct, for example, owing to a lack of available insect pollinators (33, 34) or due to increased f lower-predation (35). Interestingly, pheno- logical responses of insects also appear to be correlated with seasonal temperature (7), suggesting that plant species that respond to temperature change may better maintain important synchronous interactions, such as those between plants and pollinators (36), or better avoid negative interactions, such as predation. Alternatively, selection on f lowering phenology may be indirect by acting on phenological traits that are correlated with f lowering time (e.g., leafing out times, germination; see refs. 37 and 38). For example, earlier snowmelt in the Rocky Moun- tains has been shown to induce early spring vegetative growth in certain species, exposing young buds and f lowers to frost damage and causing declines in the sizes of some populations (39). Last, the decline of more northerly distributed species suggests yet another impact of climate change: shifting species ranges. However, in our study species range was not phylogenetically conserved, meaning that it cannot explain the phylogenetic pattern of species loss. Thus, our results suggest that f lowering- time response, and not species range, better explain the phylo- genetic nature of extinction risk among f lowering plants expe- riencing rapid climate change in Concord. For this reason, species range models that attempt to predict species response to climate change may be improved if they include species phenol- ogy, particularly the ability of species to track seasonal changes in climate. Climate change appears to have had a dramatic role in shaping the contemporary composition of the Concord f lora. Given that climate models predict at least a 1.1– 6.4 °C increase in temper- ature during this century (22), changes in the Concord f lora will likely continue to be shaped in a phylogenetically biased manner. Although phylogenetic selectivity of extinction risk has been documented in animals (22) and plants (28), our study provides the strongest evidence to date that the phylogenetic pattern of extinction risk may be due to climate change. To the extent that local extinction of species underlies their global extinction (18, 40), these results represent a link between the impacts of climate change on local community composition and broader patterns of taxonomic selectivity observed in the fossil record during past mass extinction events (41, 42). Patterns of recent species loss under rapid global climate change can potentially illuminate the processes underlying past extinction events where the pattern of loss may be well characterized, but the process is less clear (e.g., the Permian–Triassic mass extinc- tion event). In the near term, this pattern of phylogenetic selectivity is likely to have an accelerated impact on the loss of species diversity: groups of closely related species are being selectively trimmed from the Tree of Life, rather than individual species being randomly pruned from its tips. Given that climate- inf luenced loss of phylodiversity has been so great in Concord, despite 60% of the area being well protected or undeveloped since the time of Thoreau, a more global approach to conser- vation prioritization is necessary to minimize future species loss. Developing global conservation strategies will necessitate in- cluding information not only on species life history, but on their evolutionary history as well (43). Materials and Methods Study Site. Concord, Massachusetts (42°27'38'' N; 71°20'54'' W), is a small township encompassing 67 km2. Although the town has undergone extensive development since the time of Thoreau, !60% of the total area has been undeveloped or remained well protected through the efforts of numerous national, state, local, and private parks, and land-trusts (20). Floral Surveys. Thoreau surveyed the Concord area for flowering times from 1851 to 1858; Hosmer surveyed the same area from 1888 to 1902; and Primack and Miller-Rushing performed the most recent survey between 2003 and 2007 (20). Thoreau and Hosmer did not generally census graminoids, wind- pollinated trees, and wind/water pollinated aquatics due to the difficulty of determining the start of flowering; Primack and Miller-Rushing also did not sample these groups. These exclusions are not likely to affect our results for the following reasons. First, the existing sampling includes the majority (!70%) of species in Concord sensu the most comprehensive flora by Eaton (44). Second, this sampling represents all major branches of the angiosperm phylogeny (Fig. S1; www.huh.harvard.edu/research/staff/davis/Fig. S1.pdf; references for composite phylogeny construction embedded therein). Third, the exclusion of predominately wind pollinated species is not likely to have an effect on the relationship between change in abundance and flowering-time response traits: climate change appears to be much more likely to affect more conspicuously flowered, insect pollinated, species included in our dataset by means of the disruption of plant-pollinator fidelity (36). Abundance Change. The abundances of species were recorded for the 1888 – 1902 (Hosmer) and 2003–2007 (Primack and Miller-Rushing) inventories. Willis et al. PNAS ! November 4, 2008 ! vol. 105 ! no. 44 ! 17031 EV O LU TI O N http://www.pnas.org/cgi/data/0806446105/DCSupplemental/Supplemental_PDF#nameddest=ST1 http://www.pnas.org/cgi/data/0806446105/DCSupplemental/Supplemental_PDF#nameddest=ST1 http://www.pnas.org/cgi/data/0806446105/DCSupplemental/Supplemental_PDF#nameddest=ST2 http://www.pnas.org/cgi/data/0806446105/DCSupplemental/Supplemental_PDF#nameddest=SF1 Records from 1888 to 1902 included the following 6 abundance categories: Very common, common, frequent, infrequent, uncommon, and rare. Abun- dance categories from 2003 to 2007 were approximated to match the 1888 – 1902 survey by using Hosmer’s journal records, and include: very common (found throughout the area), common (occurring in (3 localities), frequent (occurring in 3 localities), infrequent (occurring in 2 localities), rare (occur in 1 locality), and very rare (10 or less individuals in a single locality). These 6 abundance categories were treated as a continuous trait scored from high (6) to low (1) abundance, with an additional scoring of zero for any species absent from a given survey. We also analyzed these data with the categories very common and common combined (i.e., states 5 to 0). This more conservative scoring did not significantly affect our results (results not shown). Change in abundance was defined as the difference in abundance between the 1888 –1902 and 2003–2007 surveys; 44 taxa that were indicated as rare in 1900 and extinct in 2007 were excluded. Rare species are considerably more likely to go extinct by chance alone (24), and so might bias our results by inflating declines in abundance. Habitat. Species were assigned to 1 of 5 habitat categories: forest, grassland and field, roadside, wetland, and aquatic. When species occurred in 2 or more habitats, they were assigned to the habitat where Eaton and Primack and Miller-Rushing saw the species most frequently (20). Habitat was included as a covariate in the models to control for the effect of habitat loss on extinction (see Phylogenetic Conservatism and Trait Correlations below). Importantly, species were lost from all habitats at approximately the same rate (20), which indicates that no habitat was particularly biased toward higher rates of extinction. This result, especially when considering the protected nature of the Concord area, indicates that these patterns of local extirpation cannot be simply explained by human development or succession. Flowering-Time Response: Flowering-Time Tracking of Seasonal Temperature. The 15-year period between 1888 and 1902 provides the longest survey period to quantify the tracking of species flowering time with seasonal temperature. Flowering-time tracking was determined with regard to seasonal variation in winter temperature (average temperature over January, April, and May; see ref. 11). April and May represent monthly temperatures commonly associated with annual flowering in this region. The month of January was also included because it was found to correlate with the flowering time of many species. This correlation is presumably due to the severe cold of midwinter, which can damage plants and, thus, delay spring flowering (23). Flowering-time tracking was quantified as the correlation coefficient between annual first flowering day and winter temperature (11). Unlike flowering-time shift, our measure of flowering-time tracking from 1888 to 1902 is less likely to be affected by changes in abundance because population size was likely more stable during this shorter period (31). This trait provides an important measure of a species ability to respond to short-term temperature variation, allowing us to relate short-term temperature response with long-term changes in abundance from 1900 to 2006. Flowering-Time Response: Shift in Flowering Time. First day of flowering was recorded by Thoreau, Hosmer, and Miller-Rushing and Primack for 465, 461, and 478 species, respectively. Observations were recorded annually for nearly all species over the duration of each botanists’ survey (11). The timing of first flowering for each species was averaged over each botanists’ survey period. Shift in first flowering day was calculated as the difference in mean first flowering day from 1850 –1900, 1850 –2007, and 1900 –2007 (11). Name Standardization. We standardized species names in the Concord flora by using the U.S. Department of Agriculture PLANTS Database (45). The most current accepted species name recognized in the database was used as our ‘‘correct’’ species name. This standardized taxonomy was then used in all downstream applications including species range estimation and phylogentic tree construction (see below). In a small number of cases (18 species), sister species were identified as synonyms. These sister taxa were collapsed into a single taxon. Species Latitudinal Range Estimation. The latitudinal data of species were compiled from several online databases including: the U.S. Department of Agriculture PLANTS Database, the National Herbarium of Canada, the Cana- dian Biodiversity Information Facility, the Royal Botanic Gardens Kew, Fair- child Tropical Botanic Garden, and the Missouri Botanical Garden (TROPICOS). Latitudinal data in these databases were derived from the literature, field- based observations, and herbarium specimens. In total, 384,292 data points were obtained for 530 species with a median of 608 observations per species. Three species with )20 observations were not included in the analysis due to the paucity of data. The average latitude for each species was obtained across the contiguous United States and adjacent Canada. The mean latitude for each species was weighted by the number of observations across the range, which more accurately represents the latitudinal affinity of each species. Local declines in species abundance could be due to populations occurring at the edge of their ranges, and thus their environmental tolerances. Alter- natively, if climate change is shifting environments northward, we would expect species with a range edge more north of Concord to be declining in abundance. We tested for the effect of species range edge on decline in abundance, and found that species with range edges north of Concord, rather than near to Concord, were much more likely to have declined in abundance. This finding supports the notion that species decline is likely associated with shifting environments resulting from climate change rather than to a local range edge effect. Because species mean latitudinal range was found to be a much better predictor of decline in abundance when analyzed with species range edge, however, the latter was excluded from our analyses. Native/Introduced Status. We obtained native/introduced status for each spe- cies from the U.S. Department of Agriculture PLANTS Database (45). Species were scored as ‘‘native’’ if they occurred in the continental United States or Canada at the time of Columbus, and ‘‘introduced’’ if they arrived from other regions since that time. A small number of species (11 species) were coded ambiguously as ‘‘native and probably introduced’’ and were not included in our analyses. Phylogeny Construction. A composite phylogeny of all species was constructed with Phylomatic version 1 (46) and was further resolved above the generic level by using recently published molecular phylogenies. Studies using (1 gene were preferred, and bootstrap support (80% was required to resolve relationships. Branch lengths were scaled to be approximately equal to time with divergence time estimates aggregated in Phylocom version 3.41 by using the ‘BLADJ’ func- tion (47). Our composite phylogeny with branch lengths scaled for time (www.huh.harvard.edu/research/staff/davis/Fig. S2.pdf, references for composite phylogeny construction embedded in Fig. S1) is available on TreeBASE (www.treebase.org). Species were pruned from this tree as necessary depending on data availability for each analysis. To test the robustness of our results to uncertainties associated with divergence time estimation, we also ran our anal- yses on the same composite tree, but with branch lengths set to 1. Phylogenetic Conservatism and Trait Correlations. The phylogenetic conserva- tism of each trait was evaluated separately by calculating the average mag- nitude of standard deviation (SD) of descendant nodes over the phylogeny, by using methods modified from Blomberg and Garland (48) as implemented in Phylocom by using the analysis of traits function (47). Standard trait correlations can be biased by species relatedness (49, 50). To account for evolutionary history in trait correlations, we used the comparative method of generalized estimating equations (GEE; ref. 51), as implemented in APE version 2.1-3 (52). GEE incorporates a phylogenetic distance matrix into the framework of a general linear model. Importantly for this study, GEE also permits the simultaneous analysis of multiple categorical and continuous traits as covariates in the same model. The inclusion of covariates allowed us to control for the effects of other factors that are likely to have an impact on change in abundance, including initial abundance, habitat, native/introduced status, and flowering season. We used 3 models to test for the correlation between change in abundance and our traits of interest (i.e., flowering-time tracking, flowering-time shift, and species latitudinal range). Model 1 tested for the effect of each trait (e.g., flowering-time tracking) on change in abundance: change in abundance & flowering-time tracking. Model 2 tested for the effect of each trait while accounting for the effects of a set of additional covariates that could also influence decline in abundance [i.e., initial abundance (ca. 1900), habitat, native/introduced status, and flow- ering season (date of first flowering)]: change in abundance ! flowering-time tracking " initial abundance " habitat " native/introduced status " flowering season. 17032 ! www.pnas.org"cgi"doi"10.1073"pnas.0806446105 Willis et al. http://www.pnas.org/cgi/data/0806446105/DCSupplemental/Supplemental_PDF#nameddest=SF1 Model 3 tested for the effect of all traits of interest (i.e., in combination) while accounting for the effects of a set of additional covariates that could also influence decline in abundance [i.e., initial abundance (ca. 1900), habitat, native/introduced status, and flowering season): change in abundance ! flowering-time tracking " flowering-time shift " species latitudinal range " initial abundance " habitat " native/introduced status " flowering season. These analyses make the assumption that intraspecific variation is less than interspecific variation. Given the phylogenetic scale at which we are compar- ing species (i.e., across all angiosperms) this is a reasonable assumption and has been demonstrated empirically (53). Sensitivity Analyses. We tested the sensitivity of our results to branch length by setting all branch lengths to 1. Also, we tested the sensitivity of our results to phylogenetic uncertainty (54). All of our analyses were tested across a set of 50 trees where the polytomies were randomly resolved on each by using the program Mesquite (55). All results were robust to these sensitivity analyses (Table S1 and Table S2). ACKNOWLEDGMENTS. We thank W. Anderson, D. Barua, A. Berry, S. Blomberg, K. Bolmgren, M. Dietze, K. Donohue, S. Edwards, K. Feeley, the Harvard Ecology Discussion Group, J. Hall, J. Hatala, R. Hellmiss, J. Losos, M. Klooster, P. Morris, B. O’Meara, L. Nikolov, S. Peters, and C. Webb for helpful discussions on this manuscript. C.C.D. and R.B.P. were supported by National Science Foundation Grants AToL EF 04-31242 and DEB 0413458, respectively. 1. Sturm M, Racine C, Tape K (2001) Climate change: Increasing shrub abundance in the Arctic. Nature 411:546 –547. 2. Beckage B, et al. (2008) A rapid upward shift of a forest ecotone during 40 years of warming in the Green Mountains of Vermont. Proc Natl Acad Sci USA 105:4197– 4202. 3. Lenoir J, Gégout JC, Marquet PA, de Ruffray P, Brisse H (2008) A significant upward shift in plant species optimum elevation during the 20th century. Science 320:1768 –1771. 4. Menzel A, et al. (2006) European phenological response to climate change matches the warming pattern. Glob Change Biol 12:1969 –1976. 5. Bertin RI (2008) Plant phenology and distribution in relation to recent climate change. J Torrey Bot Soc 135:126 –146. 6. Cleland EE, Chuine I, Menzel A, Mooney HA, Schwartz MD (2007) Shifting phenology in response to global change. Trends Ecol Evol 22:357–365. 7. Parmesan C (2007) Influences of species, latitudes and methodologies on estimates of phenological response to global warming. Glob Change Biol 13:1860 –1872. 8. Kurz WA, et al. (2008) Mountain pine beetle and forest carbon feedback to climate change. Nature 452:987–990. 9. Pounds JA, et al. (2006) Widespread amphibian extinctions from epidemic disease driven by global warming. Nature 439:161–167. 10. Fitter AH, Fitter RSR (2002) Rapid changes in flowering time in British plants. Science 296:1689 –1691. 11. Miller-Rushing AJ, Primack RB (2008) Global warming and flowering times in Thoreau’s concord: A community perspective. Ecology 89:332–341. 12. Ackerly DD (2004) Adaptation, niche conservatism, and convergence: Comparative studies of leaf evolution in the California chaparral. Am Nat 163:654 – 671. 13. Cavender-Bares J, Keen A, Miles B (2006) Phylogenetic structure of Floridian plant communities depends on taxonomic and spatial scale. Ecology 87:S109 –S122. 14. Kang H, Jang J (2004) Flowering patterns among angiosperm species in Korea: Diversity and constraints. J Plant Biol 47:348 –355. 15. Lord J, Westoby M, Leishman M (1995) Seed size and phylogeny in six temperate floras: Constraints, niche conservatism, and adaptation. Am Nat 146:349 –364. 16. Prinzing A, Durka W, Klotz S, Brandl R (2001) The niche of higher plants: Evidence for phylogenetic conservatism. Proc R Soc London Ser B 268:2383–2389. 17. Wright SJ, Calderon O (1995) Phylogenetic patterns among tropical flowering phe- nologies. J Ecol 83:937–948. 18. McKinney ML (1997) Extinction vulnerability and selectivity: Combining ecological and paleontological views. Annu Rev Ecol Syst 28:495–516. 19. Thoreau HD (1897) Walden (Houghton, Mifflin, and Company, Boston). 20. Primack RB, Miller-Rushing AJ, Dharaneeswaran K, Species lost from Thoreau’s Con- cord. Biol Conserv, in press. 21. Visser ME (2008) Keeping up with a warming world; assessing the rate of adaptation to climate change. Proc R Soc London Ser B 275:649 – 659. 22. Purvis A, Agapow PM, Gittleman JL, Mace GM (2004) Nonrandom extinction and the loss of evolutionary history. Science 288:328 –330. 23. Miller-Rushing AJ, Primack RB (2008) Effects of winter temperatures on two birch (Betula) species. Tree Physiol 28:659 – 664. 24. Lawton JH, Daily G, Newton I (1994) Population-dynamic principles. Philos Trans R Soc London Ser B 344:61– 68. 25. Lavergne S, Molina J, Debussche M (2006) Fingerprints of environmental change on the rare Mediterranean flora: A 115-year study. Glob Change Biol 12:1466 –1478. 26. Duncan RP, Young JR (2000) Determinants of plant extinction and rarity 145 years after European settlement of Auckland, New Zealand. Ecology 81:3048 –3061. 27. Bertin RI (2002) Losses of native plant species from Worcester, Massachusetts. Rhodora 104:325–349. 28. Sjöström A, Gross CL (2006) Life-history characters and phylogeny are correlated with extinction risk in the Australian angiosperms. J Biogeogr 33:271–290. 29. Charmantier A, et al. (2008) Adaptive phenotypic plasticity in response to climate change in a wild bird population. Science 320:800 – 803. 30. Davis MB, Shaw RG, Etterson JR (2005) Evolutionary responses to changing climate. Ecology 86:1704 –1714. 31. Miller-Rushing AJ, Inouye DW, Primack RB, How well do first flowering dates measure plant responses to climate change? J Ecol, in press. 32. Stenseth NC, Mysterud A (2002) Climate, changing phenology, and other life history and traits: Nonlinearity and match-mismatch to the environment. Proc Natl Acad Sci USA 99:13379 –13381. 33. Maad J (2000) Phenotypic selection in hawkmoth-pollinated Platanthera bifolia: Tar- gets and fitness surfaces. Evolution 54:112–123. 34. Kudo G, Suzuki S (2002) Relationships between flowering phenology and fruit-set of dwarf shrubs in alpine fell fields in northern Japan: A comparison with a subarctic heathland in northern Sweden. Arct Antarct Alp Res 34:185–190. 35. Parachnowitsch AL, Caruso CM (2008) Predispersal seed herbivores, not pollinators, exert selection on floral traits via female fitness. Ecology 89:1802–1810. 36. Memmott J, Craze PG, Nickolas MW, Price MV (2007) Global warming and the disrup- tion of plant–pollinator interactions. Ecol Lett 10:710 –717. 37. Primack RB (1987) Relationships among flowers, fruits, and seeds. Annu Rev Ecol Syst 18:409 – 430. 38. Donohue K (2002) Germination timing influences natural selection on life-history characters in Arabidopsis thaliana. Ecology 83:1006 –1016. 39. Inouye DW, McGuire AD (1991) Effects of snowpack on timing and abundance of flowering in Delphinium nelsonii (Ranunculaceae): Implications for climate change. Am J Bot 78:997–1001. 40. Freville H, McConway K, Dodd M, Silvertown J (2007) Prediction of extinction in plants: Interaction of extrinsic threats and life history traits. Ecology 88:2662–2672. 41. McElwain JC, Punyasena SW (2007) Mass extinction events and the plant fossil record. Trends Ecol Evol 22:548 –557. 42. Tannner LH, Lucas SG, Chapman MG (2004) Assessing the record and causes of Late Triassic extinctions. Earth-Sci Rev 65:103–139. 43. Edwards EJ, Still CJ, Donoghue MJ (2007) The relevance of phylogeny to studies of global change. Trends Ecol Evol 22:243–249. 44. Eaton RJ (1974) A Flora of Concord: An Account of the Flowering Plants, Ferns, and Fern-Allies Known to Have Occurred Without Cultivation in Concord, Massachusetts from Thoreau’s Time to the Present Day (Museum of Comparative Zoology, Harvard Univ Press, Cambridge, MA). 45. U.S. Department of Agriculture Natural Resources Conservation Service (2006) PLANTS Database (Dept. of Agriculture, Washington, DC), available at http://plants.usda.gov/. 46. Webb CO, Donoghue MJ (2005) Phylomatic: Tree assembly for applied phylogenetics. Mol Ecol Notes 5:181–183. 47. Webb CO, Ackerly DD, Kembel SW (2006) Phylocom: Software for the Analysis of Community Phylogenetic Structure and Trait Evolution. Version 3.40. Available at www.phylodiversity.net/phylocom/. 48. Blomberg SP, Garland T, Ives AR (2003) Testing for phylogenetic signal in comparative data: Behavioral traits are more labile. Evolution 57:717–745. 49. Felsenstein J (1985) Phylogenies and the comparative method. Am Nat 125:1–15. 50. Harvey PH, Pagel MD (1991) The Comparative Method in Evolutionary Biology (Oxford Univ Press, Oxford). 51. Paradis E, Claude J (2002) Analysis of comparative data using generalized estimating equations. J Theor Biol 218:175–185. 52. Bolker B, et al. (2007) APE: Analyses of Phylogenetics and Evolution. Version 2.0-2. Available at http://ape.mpl.ird.fr/. 53. Reich PB, et al. (1999) Generality of leaf trait relationships: A test across six biomes. Ecology 80:1955–1969. 54. Donoghue MJ, Ackerly DD (1996) Phylogenetic uncertainties and sensitivity analyses in comparative biology. Philos Trans R Soc London Ser B 351:1241–1249. 55. Maddison WP, Maddison DR (2006) Mesquite: A Modular System for Evolutionary Analysis. Version 1.11. Available at http://mesquiteproject.org. Willis et al. PNAS ! November 4, 2008 ! vol. 105 ! no. 44 ! 17033 EV O LU TI O N http://www.pnas.org/cgi/data/0806446105/DCSupplemental/Supplemental_PDF#nameddest=ST1 http://www.pnas.org/cgi/data/0806446105/DCSupplemental/Supplemental_PDF#nameddest=ST2 Supporting Information Willis et al. 10.1073/pnas.0806446105 Fig. S1. Composite phylogeny of 429 flowering plant species from the Concord (Massachusetts) flora depicting changes in abundance from 1900 to 2007. Branch color illustrates parsimony character state reconstruction of change in abundance as summarized in Fig. 1. For character state scoring see color legend. Lineages that exhibited an average decline in abundance of 50% or greater are indicated by a black dot. Willis et al. www.pnas.org/cgi/content/short/0806446105 1 of 3 http://www.pnas.org/cgi/content/short/0806446105 Table S1. Statistical tests of phylogenetic conservatism and trait correlations with change in abundance and with branch lengths on composite phylogeny set to 1 Trait Trait correlation Model 1 Model 2 Model 3 n Estimate n Estimate n Estimate Flowering-time tracking of seasonal temperature 175 !0.92 *** 166 !1.00 *** 140 !1.33 *** Shift in flowering-time 1850–1900 319 !0.01 * 311 !0.01 * 140 0.01 — Shift in flowering time 1900–2006 303 0.04 *** 296 0.02 *** 140 0.03 *** Shift in flowering time 1850–2006 271 0.03 *** 253 0.02 *** 140 — — Mean latitudinal range 414 !0.05 ** 362 !0.05 *** 140 !0.09 ** Change in abundance 1900–2006 — — — — — — — — — The significance of phylogenetic conservatism was tested by comparing the rank of the observed standard deviation (SD) of descendent trait means with a null model based on 9,999 random iterations of trait distributions across the composite phylogeny. The observed rank is compared with a 2-tail test of significance, i.e., an observed (obs.) rank of 250 equals a P value of 0.05. Trait correlations were tested by using a general estimator equation (GEE). Estimates describe the direction and magnitude of the correlation (e.g., a negative estimate "!0.1# of mean latitude with change in abundance suggests that species from more southerly latitudes are increasing in abundance). Model 1 (univariate model), correlation of change in abundance with each trait; Model 2 (multivariate model), correlation of change in abundance with each trait and habitat, abundance (ca. 1900), flowering season, and native/introduced status as covariates; Model 3 (multivariate model), correlation of change in abundance with all traits and habitat, abundance (ca. 1900), flowering season, and native/introduced status as covariates (shift in flowering-time response 1850 –2006 was excluded due to its high correlation with the other flowering-time shift traits). †, P $ 0.1; *, P $ 0.05; **, P $ 0.01; ***, P $ 0.001; n $ sample size. Willis et al. www.pnas.org/cgi/content/short/0806446105 2 of 3 http://www.pnas.org/cgi/content/short/0806446105 Table S2. Sensitivity analyses of phylogenetic uncertainty Trait Trait correlation Model 1 Model 2 Model 3 n Estimate n Estimate n Estimate Flowering-time tracking of seasonal temperature 175 !0.03 — 166 !1.24 *** 137 !1.58 *** Shift in flowering time 1850–1900 319 0.02 *** 311 0.01 * 137 0.01 — Shift in flowering time 1900–2006 303 0.03 *** 296 0.01 *** 137 0.01 *** Shift in flowering time 1850–2006 271 0.03 *** 253 0.01 *** 137 — — Mean latitudinal range 414 !0.01 * 362 !0.03 * 137 !0.03 — Change in abundance 1900–2006 — — — — — — — — — Phylogenetic conservatism and trait correlations tested over 50 trees with randomly resolved polytomies. Median statistic of these analyses are reported in the table. The significance of phylogenetic conservatism was tested by comparing the rank of the observed SD of descendent trait means with a null model based on 9,999 random iterations of trait distributions across the composite phylogeny. The observed rank is compared with a 2-tail test of significance, i.e., an observed rank of 250 equals a P value of 0.05. Trait correlations were tested by using a GEE. Estimates describe the direction and magnitude of the correlation (e.g., a negative estimate "!0.1# of mean latitude with change in abundance suggests that species from more southerly latitudes are increasing in abundance). Model 1 (univariate model), correlation of change in abundance with each trait; Model 2 (multivariate model), correlation of change in abundance with each trait and habitat, abundance (ca. 1900), flowering season, and native/introduced status as covariates; Model 3 (multivariate model), correlation of change in abundance with all traits and habitat, abundance (ca. 1900), flowering season, and native/introduced status as covariates (shift in flowering-time response 1850 –2006 was excluded due to its high correlation with the other flowering-time shift traits). †, P $ 0.1; *, P $ 0.05; **, P $ 0.01; ***, P $ 0.001; n $ sample size. Willis et al. www.pnas.org/cgi/content/short/0806446105 3 of 3 http://www.pnas.org/cgi/content/short/0806446105