Title: Assessing the impacts of climate change on reproductive phenology in tropical rainforests of Southeast Asia

: In humid forests in Southeast Asia, many species from dozens of plant families flower gregariously and fruit synchronously at irregular multi-year intervals 1–4 . Little is known about how climate change will impact these community-wide mass reproductive events. Here, we perform a comprehensive analysis of reproductive phenology and its environmental 5 drivers based on a monthly reproductive phenology record from 210 species in 41 families in peninsular Malaysia. We find that the proportion of flowering and fruiting species decreased from 1976 to 2010. Using a phenology model with inputs obtained from general circulation models, we show that low-temperature flowering cues became less available during the monitoring period and will further decrease in the future, leading to decreased 10 flowering opportunities in 57% of species in the Dipterocarpaceae family. Our results highlight the vulnerability of and variability in phenological responses across species in tropical ecosystems that differ from temperate and boreal biomes. biological springs and delayed biological winter arrivals are documented in response to increased temperatures 34,35 . Tropical species have been suggested to be more sensitive to climatic changes than temperate and boreal species because they have evolved in areas with less seasonal environmental variation 36,37 . However, no studies have addressed the vulnerability of tropical 30 species under future climate change. Our projections of 21st century changes in the flowering phenology of tropical plant species revealed that a 1.2°C increase in temperature under the RCP2.6 scenario resulted in an approximately 50% decrease in the future flowering probabilities of 57% of dipterocarp species that are sensitive to low-temperature flowering cues. In a temperate perennial herb that requires winter cold for floral initiation, a significantly decreased flowering 35 probability was predicted when the temperature increased by 4.5°C 38 . These results suggest that tropical species are more sensitive and vulnerable to climate change than species located in temperate ecosystems. Our results also highlight the variable features of phenological changes among species in response to climate change. Forty-three percent of dipterocarp species are predicted to be sensitive only to drought for flowering, and their reproductive activities are robust 40 to climate change. Different phenological responses across species would alter forest regeneration and, eventually, the plant species composition in the future. Regardless of the significant effect of climate change on the flowering probability at the quantitative level, the seasonal distribution of flowering probability were found to be robust in wide regions of Southeast Asia; this result represents another interesting difference between tropical and temperate plant species. climate input data from January 1, 2050, to December 31, 2099, with a daily temporal resolution and a 0.5° spatial resolution, provided by the ISI-MIP project 50 ; these data are based on the Coupled Model Intercomparison Project Phase 5 (CMIP5) outputs from three GCMs: GFDL–ESM2M, IPSL–CM5A-LR, and MIROC5. To 35 compare the flowering phenology between 1976–1996 and 2050–2099, bias-corrected GCM data from May 1, 1976, to March 31, 1996, were also used. This period (May 1, 1976–March 31, 1996) is consistent with the period used for model fitting. We selected daily minimum temperature and precipitation time series from the 0.5° grid cells corresponding to the study site for phenology monitoring at FRIM. To compare flowering phenology among regions, we also used the same set 40 of data from three other regions in Southeast Asia: Trang Province in Thailand (7° 4' N, 99° 47' E), Lambir Hills National Park in Malaysia (4° 2' N, 113° 50' E), and central Kalimantan in Indonesia (0° 06' S, 114° 0' E). Because the study site in FRIM was not in the centre of a 0.5° grid cell, we interpolated the data using four grid cells in the vicinity of the observation site. We used the weighted average according to the distance between each observation site and the centre of 45 each corresponding grid cell.

the most abundant family (45% of total species), followed by Malvaceae (7.6%) and Leguminosae (6.2%) (Fig. 1a). These long-term phenology data from more than 200 species exposed to the same environment at the arboretum provide an excellent opportunity to compare phenological responses to climate change among species.
The fractions of flowering and fruiting species fluctuated heavily between years. The 5 greatest numbers of flowering and fruiting events occurred in 1985, in which more than 35% of monitored species participated in flowering and fruiting (Fig. 1b, c). Large flowering events with the flowering of more than 20% of monitored species occurred six times over 35 studied years, and these flowering events were followed by mass fruiting events (Fig. 1b, c, Extended Data Fig.  1a). These six large reproductive events at FRIM were synchronized with general flowering events 10 monitored in natural forests in Peninsular Malaysia 8,9,23,24 , and the comparison suggests that the flowering and fruiting patterns between the arboretum at FRIM and natural forests were similar. The levels of between-species synchrony in flowering and fruiting events were significantly higher in Dipterocarpaceae species than in non-Dipterocarpaceae species (P < 0.0001, two-way analysis of variance (ANOVA), n = 4,465-6,555; Extended Data Fig. 1b). Moreover, the coefficients of 15 variation in the proportion of flowering and fruiting species in Dipterocarpaceae were twice as large as those in non-Dipterocarpaceae species (1.787 and 1.583 for Dipterocarpaceae and 0.803 and 0.753 for non-Dipterocarpaceae, respectively). These results indicate that Dipterocarpaceae plays a pivotal role in general flowering. We detected decreasing trends in the proportions of flowering (P = 0.0021, Mann-Kendall 20 (MK) test, two-sided, n = 400; Extended Data Table 1) and fruiting species (P < 0.0001, MK test, n = 400; Extended Data Table 1) from the mid-1970s to the early 2000s. In contrast, temperature and precipitation during this period revealed an increasing trend at mean rates of 0.39 ± 0.02°C (P < 0.0001, MK test, two-sided, n = 10,330; Fig. 2d) and 0.51 ± 0.26 mm/day per decade (P = 0.027, MK test, two-sided, n = 12,668; Fig. 2e), respectively. The monthly frequency of flowering and 25 fruiting varied largely among species (Fig. 2a, b); 17% of species flowered at least once per year, whereas 25% of species flowered only once every 10 years (Extended Data 2). Regardless of this large variation in the frequency of reproductive events across species, most species exhibited clear reproductive seasonality. At the community level, two flowering peaks occurred in April and October (Fig. 2a), followed by fruiting peaks in June and December, respectively (Fig. 2b). The 30 seed dispersal time, which was calculated as the month when all mature fruits dropped from their mother tree, peaked in February or August, two months after the fruiting peaks (Fig. 2c). The timing of the two seed dispersal peaks matched the phases in which temperatures and rainfall started to increase (Fig. 2d, e), which was consistent with the finding that seed dispersal is timed to match the favourable humidity condition for the survival of offspring 29 . Among the nine families 35 containing at least five species, only Moraceae, which includes the genus Ficus, produced flowers and fruits almost year-round without any seasonality (Fig. 3). Other families show synchronized flowering phenology in both spring and autumn (e.g., Dipterocarpaceae; Fig. 3) or spring flowering dominance (e.g., Anacardiaceae, Lauraceae, and Meliaceae; Fig. 3). The decreased proportions of flowering and fruiting species observed in the past can be 40 explained by increased temperatures or decreased drought event frequencies because cool temperatures and drought have been suggested to be major environmental drivers of general flowering [20][21][22] . To investigate the relationships between flowering and temperature and between flowering and precipitation, we adopted a model that was developed to predict flowering phenology in Dipterocarpaceae 21 and was further extended to predict the flowering phenology of 45 tropical plants on Barro Colorado Island, Panama 30 . The previous model was successful in explaining the weekly flowering phenology of five dipterocarp species over 10 years in the Pasoh forest reserve in peninsular Malaysia; thus, it is a reliable model for forecasting future flowering phenology in Dipterocarpaceae. The model assumed that potential environmental cues to floral induction, the cool unit (CU) and drought unit (DU), accumulate for n1 days prior to the onset of floral induction, and flowers then develop for n2 days before opening (Extended Data Fig. 2). The 5 DU model evaluates whether drought alone cues flowering (drought-induced flowering). The CU × DU model evaluates whether cold and drought have a multiplicative effect on flowering (low temperature-and drought-induced flowering). Logistic regression was performed using only DU and using CU × DU as the explanatory variables and using the presence or absence of a first flowering event as the dependent variable. We examined these two models because dipterocarp 10 flowering has been predicted successfully by both of these models in previous studies 21,29 .
To perform model fitting, we focused only on Dipterocarpaceae because of the low percentage of missing values (4.81%) compared to those of non-Dipterocarpaceae species (16.8% on average; Extended Data 1). Because some species have similar flowering phenology, we first performed time-series clustering using 98 dipterocarp species based on the similarity of their 15 flowering phenology (Fig. 4a) and then carried out forward selection of the optimal number of phenological clusters based on the minimization of the Akaike information criterion (AIC) (Extended Data Fig. 3). We used data from May 1976 to March 1996 to train the model and data from July 1997 to April 2005 to validate the model. We chose these periods for the model training and validation because there was a blank period in the dataset from April 1996 to June 1997 in 20 which climate data were missing. The model fitting results revealed that the optimal number of phenological clusters was 10 (seven clusters and three independent species; Extended Data 3). After removing independent species and clusters with fewer than 5 species (due to their small sample sizes), six clusters remained ( Fig. 4a; Extended Data 4). The CU × DU model was selected to explain the flowering 25 phenology of clusters 3 and 4, two major clusters including 27 and 28 species, respectively ( Fig.  4b; Extended Data Table 2). The flowering phenology of other clusters was explained by drought only ( Fig. 4b; Extended Data Table 2). The area under the ROC (receiver operating characteristic) curve (AUC) values ranged from 0.64-0.78 for the training data and 0.62-0.79 for the validation data (Extended Data during the 2050-2099 period across the three models decreased in clusters 3 and 4 to 57% and 49% of the predicted flowering probabilities during the 1976-1996 period, respectively (Fig. 5c). Under the RCP8.5 scenario, the mean predicted flowering probabilities in clusters 3 and 4 were further reduced to 37% and 28% during the 2050-2099 period, respectively (Fig. 5c). The decreased flowering probabilities in these two phenological clusters were caused by lower 45 occurrences of low-temperature flowering cues in the future (Extended Data Fig. 5). Under the . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 27, 2021. ; https://doi.org/10.1101/2021.08.24.457576 doi: bioRxiv preprint RCP8.5 scenario, low-temperature signals that trigger flowering in the species included in phenological clusters 3 and 4 rarely occurred or completely disappeared (Extended Data Fig. 5). In contrast, in species included in other clusters that were sensitive only to drought cues for flowering, the mean flowering probabilities were unchanged between the 1976-1996 and 2050-2099 periods (Fig. 5c) because the drought-related flowering cues were rather stable throughout 5 the simulations (Extended Data Fig. 5).
To confirm our predictions, we extended our phenology forecasting to three other regions in Southeast Asia, Trang Province in Thailand 32 , Lambir Hills National Park in Malaysia 12 , and central Kalimantan in Indonesia 13 , where long-term phenology monitoring plots exist (Fig. 6a). Decreased flowering probabilities were predicted only in clusters 3 and 4 in all regions, while the 10 flowering probabilities of other species were predicted to be robust (Extended Data Fig. 6), confirming the predictions obtained in FRIM. A comparison of seasonal flowering patterns along a latitudinal gradient based on historical climate data simulated during the 1976-1996 period from GCMs revealed shifts from a unimodal flowering peak in spring (March in the northern hemisphere) in Trang Province, bimodal or weak flowering peaks in spring and fall in FRIM and 15 Lambir National Park and a pronounced flowering peak only in spring (September in the south hemisphere) in central Kalimantan (Fig. 6b). The predicted seasonal flowering patterns are consistent with previous observations 12,13,32 , showing that using low temperatures and drought to forecast phenology can be applicable to wide regions in Southeast Asia. The latitudinal gradient of seasonal flowering patterns was predicted to be robust to climate change in the 21st century 20 (Fig. 6c), suggesting that the different seasonal distribution of flowering probability among the four regions can be explained by the differential seasonal rainfall patterns across the regions (Extended Data Fig. 7). These results suggest that the phenological responses of tropical trees in Southeast Asia to climate change are not qualitative but are rather quantitative. Phenological shifts are among the most widely studied biotic responses to climate change 33 . 25 Most phenological shift observations come from temperate and boreal biomes where advancing biological springs and delayed biological winter arrivals are documented in response to increased temperatures 34,35 . Tropical species have been suggested to be more sensitive to climatic changes than temperate and boreal species because they have evolved in areas with less seasonal environmental variation 36,37 . However, no studies have addressed the vulnerability of tropical 30 species under future climate change. Our projections of 21st century changes in the flowering phenology of tropical plant species revealed that a 1.2°C increase in temperature under the RCP2.6 scenario resulted in an approximately 50% decrease in the future flowering probabilities of 57% of dipterocarp species that are sensitive to low-temperature flowering cues. In a temperate perennial herb that requires winter cold for floral initiation, a significantly decreased flowering 35 probability was predicted when the temperature increased by 4.5°C 38 . These results suggest that tropical species are more sensitive and vulnerable to climate change than species located in temperate ecosystems. Our results also highlight the variable features of phenological changes among species in response to climate change. Forty-three percent of dipterocarp species are predicted to be sensitive only to drought for flowering, and their reproductive activities are robust 40 to climate change. Different phenological responses across species would alter forest regeneration and, eventually, the plant species composition in the future. Regardless of the significant effect of climate change on the flowering probability at the quantitative level, the seasonal distribution of flowering probability were found to be robust in wide regions of Southeast Asia; this result represents another interesting difference between tropical and temperate plant species. 45 . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made Although the data presented here are one of the longest records of reproductive phenology in tropical ecosystems 39 , we still need to be cautious when interpreting these results because there is room to extend our analyses from the phenological cluster level to the species level when longerterm and higher-temporal-resolution data become available. With upgraded phenological data, the estimation accuracy of the environmental drivers of tropical phenology and predictive power will 5 be improved. Because observations of reproductive phenology in tropical plants are still rare due to a paucity of long-term studies 26-28 , continued phenology monitoring is necessary 39 .
There have been only a few reports on past phenological changes in tropical plants. In Kibale National Park, Uganda, fruiting phenology from two data sets (covering 1970-1983 and 1990-2002) revealed that the proportion of individual fruiting was negatively correlated with the 10 minimum temperature and that increased rainfall was associated with the complete absence of fruiting in common tree species 40 . In Ranomafana National Park, Madagascar, 12 years of fruiting phenology observations showed a correlation between increased rainfall and an increase in fruiting in some tree species 41 . These studies highlight the complex and variable features of phenological changes among tropical plant species in response to climate change 37 . An improved mechanistic 15 understanding of the environmental drivers of reproductive phenology in diverse species in different tropical ecosystems will unravel the complex nature of phenological responses in the tropics and will allow the extension of future reproductive phenology projections from regional to global scales 42 .
The rapid global warming that occurred over the last millennium was unprecedented 43 . Our 20 results suggest that plastic responses to climate change at the individual level may not be high for the tropical species studied herein. Moreover, species with long generation times are the least likely to be rescued by evolution alone 44 . Early detections of biotic change signatures and predictions of the magnitude and direction of changes in plant reproductive phenology will benefit management programmes and aid in the sustainable future of the most diverse ecosystems on Earth.  . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  given for each plot.
. CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  and fruiting (grey) frequencies were drawn for nine families in which at least five species were included. The number of species included for each family is given in Fig. 1a.
. CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made   Table 3. 15 . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  1) Percentage of missing values is less than or equal to 50%: If the monthly flowering or fruiting phenology data of a species included a substantially large number of missing values (more than 25 50%), the species was excluded. 2) Stable flowering period: We considered an observation to be unreliable if the flowering period was significantly different among flowering events (if the coefficient of variation in the flowering period was larger or equal to 1.0). 3) Flowering period is shorter than or equal to 12 months: We considered an observation to be 30 unreliable if the flowering period was longer than 12 months because it was unlikely that the same tree would flower continuously for longer than 1 year. 4) The flowering and fruiting frequencies were not significantly different between the first and second half of the census period: When the flowering frequency was zero for the first half of the observation period but was larger than 0.1 for the second half of the observation period, it 35 was likely that the specimen was newly planted after monitoring began. Similarly, when the flowering frequency was zero for the second half of the observation period but was larger than 0.1 for the first half of the observation period, it was likely that the species died during monitoring. We removed these species for our analyses. We adopted the same criteria for the fruiting phenology data. 40 5) We removed overlapping species, herb species, and specimens with unknown species names.
After removing unreliable species based on the five criteria explained above, we obtained 95 dipterocarp and 115 non-dipterocarp species (Extended Data 1). We used these species for further analyses. 45 . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 27, 2021. ; https://doi.org/10.1101/2021.08.24.457576 doi: bioRxiv preprint

Detection of seasonality in reproductive phenology
To compare the flowering and fruiting phenology seasonality among different families, nine families that included at least five species were used. The number of flowering or fruiting events was counted for each month from January to December during a census, and then the frequency distribution was drawn as a histogram. Similarly, we also generated a histogram for the seed 5 dispersal month, which was calculated as the month when fruiting ended (i.e., when the binary fruiting phenology data changed from one to zero).

Classification of phenological patterns
To classify the phenological patterns, we performed time-series clustering using the R package was estimated based on AIC, as explained below. spanned a period shorter than 3 days, we approximated these missing values using the mean minimum temperatures recorded on the adjacent three days. Although solar radiation data were not available for our study, the use of precipitation is sufficient for model fitting because there is a significant negative correlation between solar radiation and precipitation in Southeast Asia 49 . 30

Climate data
Climate data generated by general circulation models (GCMs) As the future climate inputs, we used bias-corrected climate input data from January 1, 2050, to December 31, 2099, with a daily temporal resolution and a 0. . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made Although the climate input data provided by ISI-MIP were already bias-corrected, we conducted additional bias correction at FRIM using a historical scenario for each GCM dataset and the observed weather data from January 1, 1976 to December 31, 2004 based on previously presented protocol 51 . We did not implement any bias correction for the frequency of dry days or precipitation intensity of wet days 51 because we only focused on the average precipitation. 5 The variances in the annual fluctuation of the monthly mean precipitation were not the same between the observation data and historical GCM runs at FRIM. For all three GCMs (GFDL-ESM2M, IPSL-CM5A-LR, and MIROC5), the variances in the yearly fluctuation output by the GCMs tended to be larger than that of the observed data at the FRIM KEPONG weather station during winter and spring. On the other hand, during summer and fall, the variances output by the 10 GCMs tended to be smaller than that of the observed data. These biases could not be corrected using the previous method 51 . Therefore, we conducted the following bias correction for these data: where #,$ %&' is the normal climate value during the target period. In this method, we defined the normal climate value as the mean of the monthly mean precipitation values over 31 years. (4) 35 When the mean of a gamma distribution is fixed at one, the shape parameters are represented as follows: 40 . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made In this method, we assumed that the #,$ %&' value follows a gamma distribution and that the ratio of the variance of #,$ %&' to the variance of #,$ 9:; is maintained even in the future scenario.
Here, #,$ 9:; represents the deviation of the monthly mean in the observation data from the normal 5 climate value. (8) 10 In the above equations, ̅ #,$ 9:; indicates the monthly mean precipitation value in the observed data.
As mentioned above, because we assume that the ratio of the variance in #,$ %&' to the variance in (11) 20 In equation 11, #,$ %&' , is the deviation of the monthly mean of the historical GCM precipitation data from the normal climate value. Here, we defined the normal climate value as the average monthly mean during 1976-2004.
The method proposed here is an original bias correction method, but the above equations 25 are easily derived if we assume that the #,$ %&' value follows a gamma distribution and that the ratio of the variance in #,$ %&' to the variance in #,$ 9:; is maintained even in the future scenario. Notably, because we combined this method with the bias correction method described previously 51 , equation 2 should be expressed as follows: 30 where 8 G,#,$ %&' is the precipitation data that are bias-corrected using the method described previously 51 . 35 Analyses We adopted previously presented 21 models in which environmental triggers for floral induction accumulate for n1 days prior to the onset of floral induction (Fig. 3A). Flowers then develop for n2 days before opening (Extended Data Fig. 2). The model assumption of the time lag between floral . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 27, 2021. ; https://doi.org/10.1101/2021.08.24.457576 doi: bioRxiv preprint induction and anthesis, which is denoted as n2, was validated by a previous finding in which the expression peaks of flowering-time genes, which are used as molecular markers of floral induction, were shown to occur at least one month before anthesis in Shorea curtisii 21 . S. curtissi is included in our data set. The cool unit at time t, CU( | H ), is calculated as follows: where H = { + , < , ̅ } is the set of parameters and x(t) is the temperature at time t. Here, ̅ indicates the threshold temperature. The term max{x1, x2} is a function that returns a larger value for the two arguments. Similarly, given J = { + , < , K }, the drought unit at time t, DU( | J ), is defined as the difference between the mean daily accumulation of rainfall over n1 days and a 10 threshold rainfall level ( ): where y(t) is the rainfall value at time t. The term max{x1, x2} is defined similarly as in equation 13.

15
Logistic regression was performed using only the drought unit (DU) and using the product of CU and DU (CU × DU) as the explanatory variables and using the presence or absence of a first flowering event as the dependent variable for each phenological cluster. Because the number of phenological clusters is unknown, we performed forward selection on the cluster number based on the AIC. Let m be the number of phenological clusters based on the dendrogram drawn from the 20 time-series clustering explained above (Extended Data Fig. 5). Given m phenological clusters, let where the element in the bracket indicates the ID of the cluster in which the DU model is adopted for model fitting. When k = 0, the DU model is not used; instead, the CU × DU model is adopted for model fitting for both clusters 1 and 2. Let i be the ith element of the vector E, which is defined 30 as follows: where n is the length of the time series data for each cluster. Notably, n = 223 is the same for all (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 27, 2021. ; https://doi.org/10.1101/2021.08.24.457576 doi: bioRxiv preprint where #,0 ( ) is the dummy variable indicating a cluster for i; #,0 ( ) equals 1 if the ith element of E belongs to the jth cluster, otherwise it is zero, and #,0 and #,0 in equation (5) are regression coefficients for the jth cluster when the species are grouped into m clusters. We estimate the parameters and the number of clusters based on a finite number of observations. Given the number of clusters m, for each of m clusters, the parameters were estimated by maximizing the 5 loglikelihood value calculated for all combinations of potential parameter values for + , < , ̅ , and K within the ranges of [1 (min), 50 (max)] for n1, [1, 50] for n2, [19,25] for ̅ , and [1, 9] for K . We varied the days (n1 and n2) by integers, temperature ( ̅ ) by tenths of a degree C, and daily precipitation ( K ) by tenths of a mm. Regression coefficients ( #,0 , #,0 ) for all j values under a given m value and associated likelihoods were determined using generalized linear models with 10 binomial error structures.
With the results of the parameter estimations, we determined the number of clusters in two steps. For the first step, for a given m, we obtained T ( ) according to the following equation: For the second step, with the results of T obtained from the first step, we obtained the estimate of the number of clusters according to forward selection by searching for the " value that satisfies the following inequalities: 20 AIC( " , T ( ")) < AIC( " + 1, T ( " + 1)) ∩ AIC( ", T ( " )) < AIC( " − 1, T ( " − 1)).
For model fitting, the first flowering month was extracted from the flowering phenology data. When flowering lasted more than one month, the month after the first flowering month was replaced by a value of zero (absence of flowering). If the month before the first flowering month 25 was a missing value, the first flowering month was treated as a missing value and was not used for further analyses. We assumed that phenology monitoring was performed on the first date of each month.

Projections of 21st century changes in flowering phenology 30
We used two scenarios (RCP2.6 and RCP8.5) to forecast future reproductive phenology in dipterocarp species for each of the three GCMs (GFDL-ESM2M, IPSL-CM5A-LR, and MIROC5). We predicted the flowering probability per month for each phenological cluster during the periods from May 1, 1976-March 31, 1996 and from January 1, 2050-December 31, 2099 based on the best model (Extended Data Table 3). The predicted flowering probability during the 35 2050-2099 period was normalized to that during the 1976-1996 period for each climate scenario and for each of three GCMs. To compare the seasonal patterns between 1976-1996 and 2050-2099, the predicted flowering probability was averaged for each month from January to December and plotted for each month in Fig. 6. R version 3.6.3 52 was used for all analyses. 40 Data availability The phenology data used in this study will be provided as the Extended Data after acceptance of the paper. The meteorological data provided by the Malaysian Meteorological Department may not be resold, distributed or disclosed to third parties, as regulated by the Malaysian . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 27, 2021. ; https://doi.org/10.1101/2021.08.24.457576 doi: bioRxiv preprint Meteorological Department. Other data that support the findings of this study are available from the corresponding authors (SN and AS) upon reasonable request.

Code availability
The codes used for model fitting will be provided as Source Data after acceptance of the paper. . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  The cool unit (CU) was equal to the cumulative sum of the negative differences between the threshold temperatures ( ̅ ) and the daily mean temperatures over n1 days. The drought unit (DU) equalled the difference between the threshold precipitation level ( K ) and the cumulative sum of precipitation over n1 days. Two models were used to describe the relationship between the 10 flowering probability and environmental cues. The two models involved the drought only (DU) and synergistic cold plus drought (CU × DU) models.
. CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  m is the number of phenological clusters examined for model fitting. The species grouped into the same cluster at m = 10 are illustrated by the same colour. The optimal cluster number ( " = 10) was identified according to the forward selection of optimal " values based on the minimization of the AIC. When m = 11, cluster 3 was divided into two subclusters (clusters 3-1 and 3-2). Clusters with fewer than four species (cluster 7 and independent species sp. 36, sp. 63 and sp. 72) were 10 removed before phenology forecasting was conducted due to their small sample size.
. CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made    Fig. 6a).
. CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made . CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

10
. CC-BY-NC-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 27, 2021. ; https://doi.org/10.1101/2021.08.24.457576 doi: bioRxiv preprint