Mechanisms of changes in spring phenology

Plant phenology, particularly in spring, is shifting with warming climate (Fitter and Fitter, 2002; Post et al., 2018; Primack and Gallinat, 2016). Numerous observational studies and controlled experiments, conducted over a range of climatic and phenological conditions (Huang et al., 2020; Parmesan, 2007; Prevéy et al., 2019; Root et al., 2003; Wolkovich et al., 2012), have produced varying results, with observed phenological changes differing by research approach (Wolkovich et al., 2012), climatic region (Parmesan, 2007; Post et al., 2018; Prevéy et al., 2017; Zhang et al., 2015), and functional group (Parmesan, 2007; Root et al., 2003; Willis et al., 2010; Wolkovich et al., 2013; Zettlemoyer et al., 2019). For example, observed changes are often smaller with experimental studies (Wolkovich et al., 2012), native species (Wolkovich et al., 2013), and warmer regions (Prevéy et al., 2017). Understanding these differences (also referred to as discrepancies, disparities, or mismatches) is crucial for predicting changes in plant communities and ecosystems (Fitter and Fitter, 2002; Polgar et al., 2014; Primack and Gallinat, 2016) but difficult using sensitivity analysis that is based on rates of phenological changes, not on drivers of spring phenology.

Plants in the northern hemisphere require an accumulation of cool winter temperatures (winter chilling) to break dormancy and an accumulation of warm spring temperatures (spring forcing, degree days) to initiate budburst (e.g., leafing or flowering) (Piao et al., 2019; Way and Montgomery, 2015). If chilling requirements are fully met and other phenological constraints (e.g., photoperiod effect and environmental stresses) do not change with warming, the changes in spring phenology, in response to a warmer climate, result primarily from changes in spring temperatures, i.e., changes in spring forcing and rate of forcing accumulation (i.e., temperature) at budburst with the warmer climate. That is, higher budburst temperatures are associated with smaller observed changes and sensitivity, opposite to the effects of forcing changes (Chu et al., 2021, 2023). As species often differ in timing of budburst, both forcing change and budburst temperature are species-specific, reflecting variations of individual species in phenological responses to climate warming (Chu et al., 2021, 2023). The effects of insufficient winter chilling, photoperiod (day length), and environmental stresses (e.g., drought and spring freezing) are similar in constraining the advance of spring phenology (Chen et al., 2011; Huang et al., 2019; Körner and Basler, 2010; Ma et al., 2019; Way and Montgomery, 2015). Therefore, the difference between observed phenological change and that expected from forcing change and budburst temperature represents phenological difference due to changes in phenological constraints and can be thereafter named phenological lag. Separating effects of different constraints is possible if soil moisture/plant water status is monitored (Huang et al., 2019; Post et al., 2022), or species-specific chilling and photoperiod effects (Man et al., 2017a, 2021a; Fu et al., 2019) or cold hardiness (Man et al., 2017b, 2021b) are known from previous research. For example, phenological lag can represent chilling effects when photoperiod effect and environmental stresses do not change with warming (Chu et al., 2021).

The changes in spring phenology are often assessed by sensitivity, a parameter that may reflect rates of phenological changes but not uneven warming effects and responses (Beaubien and Hamann, 2011; Keenan et al., 2020; Post et al., 2018; Rafferty et al., 2020), nor does the parameter provide meaningful indications of magnitude and underlying mechanisms of changes in spring phenology or differences in phenological responses (Chu et al., 2021, 2023). For example, the same average temperature change (spring or annual temperatures) with uneven warming can have greater (greater warming prior to budburst) or smaller (smaller warming prior to budburst) phenological change than that with even warming. Under these circumstances, interpretations based on sensitivity analysis are likely inappropriate. Here, we present a new method to partition observed phenological changes based on changes in spring phenology and temperatures with control (baseline) and warmer climates. We applied this analytical approach to meta-analyze phenological changes reported in the literature and investigated the mechanisms of the differential responses reported previously on research approach (observations or experiments), species origin (native or exotic), climatic region (boreal or temperate), and growth form (tree, shrub, herb, or grass). Our objective was to determine how phenological responses differ among different groups and how differential responses are related to the drivers of spring phenology, i.e., forcing change, budburst temperature, and phenological lag.

Partitioning observed changes

Species-specific forcing change (FC) can be calculated from the difference in degree days (> 0 °C, see Man and Lu, 2010) between baseline (or control) and warmer climates at budburst (leafing or flowering) with the baseline climate (OCi) (see Figure 1).

A diagram showing the relationships among forcing change (FC, difference in spring forcing (degree days) between baseline and warmer climates at budburst with the baseline climate OCi), expected response (NE, difference between baseline and warmer climates in reaching species forcing needs with the baseline climate, i.e., spring forcing at OCi), budburst temperature (FC/NE, average temperature or rate of forcing accumulation at budburst with the warmer climate), and phenological lag (NC, difference between expected NE and observed NO responses with the warmer climate).

The phenological response expected under the null hypothesis that climate warming does not induce changes in phenological constraints can be determined from forcing change and species phenology with the baseline climate (Figure 1):

where NE is the difference in the number of days (n) between baseline and warmer climates in reaching species forcing needs with the baseline climate (spring forcing at OCi). Under this hypothesis, spring phenology is expected to shift from OCi with the baseline climate to EWi with the warmer climate. NE is typically positive, meaning earlier leafing and flowering, but can be negative if temperatures decrease over time in observational studies.

Budburst temperature (TB) is the average temperature or rate of forcing accumulation within the window of expected phenological response (between OCi and EWi) and can be determined from forcing change and expected response.

The difference between expected (NE) and observed (NO) responses is phenological lag (NC) reflecting phenological difference due to change in phenological constraints induced by warming.

NC > 0 indicates increasing phenological constraints that may result from increasing chilling deficiency, photoperiod restriction, moisture stress, or risks of spring frosts; NC = 0 indicates no changes; and NC < 0 indicates decreasing constraints with warmer climate.

Meta-analysis of global data

Phenological Data

We followed the guidelines of PRISMA (preferred reporting items for systematic reviews and meta-analyses) (O’Dea et al., 2021; Page et al., 2021) to build up the dataset for responses of plant spring phenology to experimental and climatic warming. We used Web of Science and Google Scholar to search experimental and observational studies on warming-induced changes in spring phenology. The following criteria were used to select papers: (a) accessible peer-reviewed articles published in scientific journals for boreal (MAT < 6 °C) and temperate (MAT ≥ 6 °C) regions in the northern hemisphere, (b) studies reporting spring budburst (leafing or flowering) before Julian day 213 (July 31), (c) studies containing different climatic conditions of baseline (lower temperatures in early period of observations or unwarmed control treatment) and warmer (higher temperatures in later period of observations or warmer treatment) climates, (d) access to local temperature data (only one study was excluded due to lack of nearby temperature data within 50 km), and (e) reported phenological changes in unique locations, species, and periods of observations/warmer treatments. When publications did not contain sufficient phenological information, we accessed data directly from phenological databases (e.g., China’s National Earth System Science Data Center). When results were reported in graphical format, we extracted the required data using the distance measuring tool in Adobe Acrobat Reader DC. In experimental studies with multiple warming treatments, we selected the treatment with the smallest temperature increases to represent a likely scenario of climate change and to maintain data independence. We also excluded studies where the timing of spring budburst, total temperature change, or observational time periods were not available.

In total, 66 studies were identified from 87 locations, containing 1,527 phenological responses for 980 species (Figure S1; Tables S1–S3). For studies reporting multiple phenological stages, we only used data for the first occurrence of the events. To account for variation in phenological responses, each study location was characterized by research approach (observation or experiment), climatic region (boreal or temperate), biogeographic origin of species (exotic or native to study area), and growth form (tree, shrub, herb, or grass).

Climate Data

Baseline and warmer temperature data for calculating forcing change, expected response, budburst temperature, phenological lag, spring warming (average temperature change from Julian day 1 to 182), and long-term MAT and MAP were obtained from 9 different sources (Table S4). When temperature data were available from multiple sources, we selected the data set with the fewest missing values. For each phenological response, the number of weather stations within 50 km of data collection ranged from 1 for studies at a single location to 14 for observational studies across large geographic areas (Table S2).

Data analysis

We calculated forcing change, expected response, budburst temperature, and phenological lag for each of the 1527 responses compiled. We used linear mixed-effects models to explore the differences in observed responses and phenological lags between research approaches (observations vs. experiments), species origins (native vs. exotic), and climatic regions (boreal vs. temperate), and among growth forms (trees, shrubs, herbs, and grasses). Observed response is the total phenological change observed (commonly called phenological change or response), while phenological lag represents the lag effect of all constraints and therefore modification of phenological change from the expectation based on changes in spring temperatures. Separate analysis was conducted for each of the four fixed effects and two budburst events (leafing and flowering). Location and species were treated as random effects.

To assess the influences of climatic, phenological, and biological variables, we used a stepwise regression with automated combined forward and backward selection by Akaike information criterion (AIC) (Burnham and Anderson, 2002) to select the best combination of variables for predicting observed responses (NO).

The variables included were altitude (AL), latitude (LA), MAT, MAP, spring phenology (DT), budburst temperature (TB), forcing change (FC), spring warming (SW), and phenological lag (NC), and not standardized prior to regression analysis.

Observed changes in spring phenology

Across 66 studies with 980 species and 1527 responses, observed responses (NO) ranged from −23.8 to +26.8 days for leafing and −36.3 to +49.9 days for flowering, with an average of +6.0 days for both events (Table S3). These changes differed from expected responses (NE) which ranged from −3 to +27 days for leafing and −3 to +26 days for flowering and averaged +7.8 and +9.5 days, respectively. Observed phenological responses were thus smaller than the expected by 1.8 and 3.5 days, on average, for leafing and flowering, respectively.

In experimental studies, observed flowering response was non-significantly smaller than those from observational studies (1.7 days, p=0.232, Figure 2a), while phenological lag was significantly longer for both leafing (3.0 days, p=0.028) and flowering (2.7 days, p=0.027) (Figure 2b). Observed flowering response in exotic species was 2 days greater than in native species (p=0.016), likely due to differences in phenological lag (1.9 days, p=0.019). Observed leafing and flowering responses in boreal region were 4.3 and 1.9 days smaller (p=0.004, p=0.088), respectively, than in temperate region, partially due to smaller forcing changes (equivalent to 2.8 days in leafing and 1.4 days in flowering based on forcing changes and budburst temperatures, Table 1 note a). Observed leafing response in grasses was 3.8 days smaller than in trees and shrubs (p=0.033), likely due to differences in forcing changes (6.0 days, Table 1 note b).

Observed responses and phenological lags (least square means ± standard errors) in leafing and flowering (left and right set of bars in each pair, respectively) by research approach (observation and experiment), species origin (native and exotic), climatic region (boreal and temperate), and growth form (tree, shrub, herb, and grass) extracted from reported plant phenological changes in spring.

Phenological lags are calculated from the differences between observed responses and those expected from species-specific changes in spring temperatures. Different letters indicate means differ significantly (p<0.05).

Mean values of climatic, phenological, and biological variables by budburst (leafing or flowering), research approach (observation and experiment), species origin (native and exotic), climatic region (boreal and temperate), and growth form (tree, shrub, herb, and grass).

Budburst temperature =average temperature at budburst with the warmer climate, forcing change =warming-induced changes in spring forcing prior to budburst, spring warming=average spring temperature change, and phenological lag= difference between observed response and that expected from forcing change and budburst temperature.

Other than phenological lag, forcing change and budburst temperature were the most influential variables affecting observed responses (Table 2). The remaining variables combined (altitude, latitude, MAT, and spring warming) explained <2.5% variations in observed leafing and flowering responses.

Final stepwise regression coefficients and variable influence on observed phenological responses by Akaike information criterion (AIC).

Differential phenological responses

Observations vs. experiments

By partitioning observed changes based on drivers of spring phenology, we clarified some of the uncertainty regarding the underlying mechanisms of differences in plant phenological response to climate warming. By synthesizing global data, we found that observed flowering responses in experimental studies were non-significantly smaller than that in observational studies, but this was not the case for leafing responses; these findings differ from those based on temperature sensitivity (Wolkovich et al., 2012). The greater warming in experimental studies did not produce greater observed responses, due to longer phenological lags, which may have resulted in the reduced sensitivity (Wolkovich et al., 2012). First, experimental studies are often conducted at higher altitude sites with lower MAT, later spring phenology, and therefore longer accumulation of winter chilling. The longer phenological lags are therefore unlikely the result of insufficient winter chilling. Second, artificial warming has typically provided higher forcing changes and budburst temperatures than natural climate warming (see Table 1). Greater forcing changes produce larger expected responses (observed response + phenological lag), whereas higher budburst temperatures reduce expected responses (Chu et al., 2021; Prevéy et al., 2017; Wolkovich et al., 2021). The high temperatures in the artificial warming structures (Marion et al., 1997) can also reduce humidity and soil moisture (Ettinger et al., 2019; Huang et al., 2019), slowing spring development (Ganjurjav et al., 2020; Huang et al., 2019; Moore et al., 2015), particularly in early spring when temperatures in the surrounding environment are low, thus restricting soil water movement. While the different forcing changes adequately accounts for the differences in expected leafing responses (2.5 out of 3.3 days, Table 1 note c), the large differences in forcing changes but similar expected responses in flowering (8.63 and 9.07 days, Figure 2) suggests a strong effect of budburst temperature (Table 1 note d) (Chu et al., 2023). Thus, warming-induced moisture stress and high budburst temperatures may interact to alter observed phenological responses and temperature sensitivity in experimental studies (Wolkovich et al., 2012). Future experiments should consider humidity, soil moisture, and plant water status to better understand disparities between observational and experimental studies and relate plant phenological changes from experimental warming to those from natural climate change (Wolkovich et al., 2012).

Native vs. exotic species

Exotic species are often reported with more noticeable observed responses to warming, a trend that is generally based on flowering (Calinger et al., 2013; Willis et al., 2010; Wolkovich et al., 2013). The ecological implications of this response are well described (Polgar et al., 2014; Zettlemoyer et al., 2019), but the underlying mechanisms are not well understood. Lower budburst temperatures faced by early-start exotics are likely partially the cause (Chu et al., 2021, 2023), as shown by a reverse trend for late-start exotics (Zohner and Renner, 2014). In this synthesis, phenological lag did not increase from early leafing to late flowering for exotic species, a trend that is consistent with dormancy release commonly reported on leaf and flower buds (Campoy et al., 2013; Gariglio et al., 2006; Hussain et al., 2015; Wall et al., 2008), contrary to that for native species (Figure 2b). The smaller observed response in flowering for native species is associated with longer phenological lags and accumulations of winter chilling (Table 1), suggesting more stressful environments later in the season or higher stress sensitivity with reproductive events. Comparatively, plants starting growth early in the growing season may be less restricted by soil moisture that often recharges from winter precipitation and is depleted with increased moisture consumption over time (Sherry et al., 2007; Wolkovich et al., 2013; Zettlemoyer et al. 2019), particularly in dry climates (Man and Greenway, 2013). Moisture stress can progressively develop over the season (Piao et al. 2019; Wolkovich et al., 2013; Zettlemoyer et al., 2019) and differentially affect exotic and native species with differing spring phenology (Stuble et al., 2021; Willis et al., 2010; Wolkovich et al., 2013) The smaller observed responses in flowering of native species suggests reproductive disadvantage (IPCC, 2014; Zettlemoyer et al., 2019). The consistent leafing response, however, does not support suggestions that late-start native species would have relatively shorter active growing seasons and, therefore, a competitive growth disadvantages relative to early-start exotic species with climate warming (Fitter and Fitter, 2002; Morin et al., 2009; Primack and Gallinat, 2016). Total thermal benefits from climate warming do not differ among species with differing spring phenology (Chu et al., 2021). Given the projected increases in temperature and decreases in precipitation for many parts of the world (IPCC, 2014), future studies should assess the differences in drought sensitivity between vegetative and reproductive events and among functional groups and climatic regions.

Boreal vs. temperate region

Spatial variations along altitudinal and latitudinal gradients are often reported but not well understood (Ge et al., 2015; Morin et al., 2009; Parmesan, 2007; Zhang et al., 2015). Our analysis indicates smaller observed responses in boreal region (MAT<6 °C) due to less forcing changes (Tables 1), but not to longer phenological lag (Figure 2). The similar spring warming but less forcing changes suggests that the temperature increases in boreal region occur more frequently in the winter (Beaubien and Hamann, 2011; Shen et al., 2015) when temperatures are below freezing and don’t contribute to spring forcing (Man and Lu, 2010). The uneven warming has been reported in both natural and experimental settings (Beaubien and Hamann, 2011; Prevéy et al., 2017; Shen et al., 2015; Slaney et al., 2007; Yang et al., 2020) and supports the claim that average temperature changes do not represent species-specific forcing changes (Chu et al., 2021; Keenan et al., 2020; Wolkovich et al., 2021). When forcing changes are held constant, however, observed responses increase with altitude and latitude, as shown by the positive regression coefficients in both leafing and flowering models (Table 2). The possible mechanisms behind this may be lower budburst temperature and less moisture stress at higher altitude and latitude (Rafferty et al., 2020), as shown by the negative collinearity of altitude and latitude with budburst temperature and phenological lag in both leafing and flowering (data not shown). Therefore, boreal region, depending on forcing changes, can have smaller (Ge et al., 2015; Menzel et al., 2006; Rafferty et al., 2020; Shen et al., 2015) or greater (Ge et al., 2015; Parmesan, 2007; Post et al., 2018; Prevéy et al., 2017) phenological responses, spatial trends that may be different from those suggested by phenological sensitivity (Liu et al., 2019; Post et al., 2018). Consequently, climate change may not result in general phenological convergence across altitudes and latitudes (Rafferty et al., 2020; Shen et al., 2015; Tao et al., 2021), contrary to suggestions by others based on small scale studies (Prevéy et al., 2017; Vitasse et al., 2018; Ziello et al., 2009).

Trees, shrubs, herbs, and grasses

Smaller observed responses in herb and grass leafing support findings from some early syntheses (König et al, 2018; Parmesan, 2007), but not those that show a greater response from non-woody plants (Ge et al., 2015; Post and Stenseth, 1999; Root et al., 2003) or no differences (Stuble et al., 2021). Herbs and grasses often resume growth earlier than woody shrubs and trees (Badeck et al., 2004; Heberling et al., 2019) and should have lower budburst temperatures and therefore greater phenological responses to climate warming (Chu et al., 2021; Prevéy et al., 2017; Wolkovich et al., 2021). However, different growth forms are often studied separately on sites in different climatic regions, making comparisons challenging. Furthermore, most studies on herb and grass leafing are conducted in boreal region at high altitudes, low MAT, late spring phenology, low forcing changes, and low spring warming (Table 1); the extremely small forcing changes resulted in small observed responses despite lower budburst temperatures (Figure 2a). The observed responses with herbs and grasses would not be smaller if converted to temperature sensitivity (see Table 1 for large differences in spring warming for leafing). In contrast to leafing, herb and grass flowering studies have greater forcing changes and higher budburst temperatures (Table 1), the latter would lower both expected and observed responses (Table 1 note d) (Chu et al., 2021; Prevéy et al., 2017; Wolkovich et al., 2021) and potentially induce moisture stress and increases phenological lags (Ganjurjav et al., 2020; Huang et al., 2019; Moore et al., 2015; Post et al., 2022), particularly for grasses (Sherry et al., 2007; Zettlemoyer et al., 2019) that have shallower root systems (Kulmatiski and Beard, 2013; Schenk and Jackson, 2002).

Mechanistic assessment of changes in spring phenology

Averaged across all observational and experimental studies, observed responses and phenological lags are both positive, suggesting that climate warming advances spring phenology but not at rates expected from changes in spring temperatures. The positive lags are likely due to more stressful environments with warmer and drier climate (Huang et al., 2019; Sherry et al., 2007; Zettlemoyer et al., 2019). We quantified phenological lags from changes in spring phenology and temperatures, an approach that does not require species-specific chilling or forcing needs that are often unavailable (Fitter and Fitter, 2002; Morin et al., 2009) or variable in methods of estimation across species and studies (Man and Lu, 2010; Ettinger et al., 2020; Zhang et al., 2018). The relatively early stage of climate warming and association of longer phenological lags with longer accumulation of winter chilling (Table 1) suggest that current climate warming is not likely to induce a general chilling shortage (Chu et al., 2023; Ettinger et al., 2020; Tao et al., 2021; Yang et al., 2020), although incidental effects on individual species or in particular conditions can occur (Fu et al., 2015). Similarly, the variable and minor effects of photoperiod (Basler and Körner, 2012; Chuine et al., 2010; Ettinger et al., 2020; Zohner et al., 2016) (insignificant spring phenology; see Table 2) often reported in studies with extreme warming scenarios (Caffarra and Donnelly, 2011; Fu et al., 2019; Laube et al., 2014) or compounded with that of budburst temperatures (Chu et al., 2021), minor effects by nutrient availability (Piao et al., 2019), or incidental spring freezing (Man et al., 2009, 2021b) are not likely to cause the systemic differences in phenological lag between observational and experimental studies or between native and exotic species.

analysis indicates that both forcing change (quantity) and budburst temperature (rate of forcing accumulation) strongly influence the magnitude of change in spring phenology, consistent with the findings by Chu et al. (2021) that smaller phenological response with late season species are largely due to their higher budburst temperatures. In phenological research, the influence of budburst temperature is not adequately recognized (Chu et al., 2021, 2023) and can be confounded with progressive increases of phenological lag with greater chilling requirements or insufficient winter chilling (Asse et al., 2018; Morin et al., 2009; Primack and Gallinat, 2016), photoperiod constraints (Körner and Basler, 2010; Shen et al., 2014; Way and Montgomery, 2015), or moisture stress (Piao et al., 2019; Wolkovich et al., 2013; Zettlemoyer et al., 2019) suggested for late season species. The difference of the significance in budburst temperatures between leafing and flowering models (Table 2) suggests that the influence of budburst temperatures can be greater if budburst temperatures decrease with the advance of spring phenology and progress of climate warming (Tao et al., 2021). Compared to spring average temperature change that explained <0.5% of the variation in observed responses (Table 2), forcing change quantifies climate warming effect relevant to phenological change of individual species and reduces uneven warming effects among species with different spring phenology (Chu et al., 2021) and areas in different climatic regions (Post et al., 2018; Yang et al., 2020). Both forcing change and budburst temperature can be readily extracted from phenological and temperature data, as demonstrated in this synthesis, but have not been used in assessment of plant phenological changes in spring.

Conclusions and caveats

In this article, we outline an analytical framework to partition phenological changes based on drivers of spring phenology and report differential phenological responses identified through the meta-analysis of observed changes and phenological lag. Longer phenological lag, likely resulting from more stressful environments with warmer and drier climate, helped to explain smaller responses with experimental studies and native plants in flowering, while less forcing changes were mainly responsible for the smaller responses in leafing and flowering in boreal region and in grass leafing. Some of these differential responses are different from those reported previously based on sensitivity analysis.

Our approach does not require species-specific chilling and forcing needs, chilling– forcing relationships, or temperature response models that are often unavailable or variable in methods of estimation across species and studies. As chilling and forcing models reflect physiological aspects of chilling and forcing processes and would not change with climate warming, the use of alternative base temperatures or forcing models would not affect the partitioning of phenological changes, except for derived budburst temperature that represents the rate of forcing accumulation not actual temperatures.

While our method helps to understand changes in spring phenology and differences in plant phenological responses across different functional groups or climatic regions, the ecological implications of phenological lag can be uncertain without investigation of individual constraints. In this synthesis, the effects of photoperiod and insufficient winter chilling are likely limited across all studies with the current level of climate change. In the boreal region with a long winter, insufficient winter chilling is unlikely to occur with the levels of climate warming projected (Chu et al., 2021; Tao et al., 2021). Phenological lag indicates the overall lag effect and the need to investigate contributions of individual constraints that can be specifically determined at individual study level if biological and environmental constrains are known from on-site monitoring or previous research.

Supporting Information

Temperature data and R codes for calculating forcing changes, expected responses, budburst temperature, spring warming, and phenological lags of individual studies are deposited at Dryad https://doi.org/10.5061/dryad.dncjsxm9x.

Acknowledgements

This work was supported by Ontario Ministry of Natural Resources and Forestry (Canada), Guangxi Normal University (China), Jilin Provincial Academy of Forestry Sciences (China), Shanghai Botanical Garden (China), Research Institute of Subtropical Forestry (China), and Lakehead University (Canada). Lisa Buse and Gillian Muir of Ontario Ministry of Natural Resources and Forestry and three anonymous reviewers provided constructive suggestions for improving earlier drafts of the manuscript. Phenological and temperature data provided by China’s National Earth System Science Data Center are greatly appreciated.

Additional information

Author Contributions

R.M. and J.T. conceived the general idea of the method; Y.J., X.C., X.Y., R.M., and J.T. compiled data; and Y.J., S.J.M., X.C., X.Y., R.M., J.T., and Q.L.D. contributed to the interpretations of the results and to the writing of the paper.

Additional files

Tables S1, S2, S3, S4

Figure S1