An ecological transcriptome approach to capture the molecular and physiological mechanisms of mass flowering in Shorea curtisii

Study species and sample collection

The species selected for this study, S. curtisii (Dipterocarpaceae), and other trees from the genus Shorea section Mutica, have been documented as reliable indicators of GF (Ashton, Givnish & Appanah, 1988). Shorea spp. can be found occupying ridges of dipterocarp forests in the Malay Peninsula ranging from 300–800 m above sea level (Appanah & Chan, 1981).

Due to the unique challenges associated with sampling from the canopy of these trees (∼40 m above the ground canopy), which requires the setup of ladders along the main trunk, we were limited to accessing and collecting samples from canopies of only two S. curtisii individuals. These trees, designated as C1 and C2, had diameters at breast height of 72 cm and 89 cm, respectively (Yeoh et al., 2017). The trees are located at Semangkok Forest Reserve (2°58′N, 102°18′E, 340–450 m above sea level) in the state of Selangor Darul Ehsan, Peninsular Malaysia (Fig. S2 A). The permission to conduct this research in the forest reserve was granted by the Selangor State Forest Department. The forest reserve consists of a 6-ha (200 m × 300 m) hill dipterocarp forest conservation plot, and another 5.4-ha (140 m × 400 m) plot established on a selectively logged area (Tani et al., 2015; Yagihashi et al., 2010). Daily temperature and precipitation data were collected from Hospital Kuala Kubu Bharu meteorological station and Kampung Pertak hydrological station, respectively (Fig. S2B), and were previously reported in Yeoh et al. (2017). The cumulative rainfall data (defined as 30-d moving rainfall total) was calculated from the daily rainfall record. Drought is defined as <40 mm cumulative rainfall (Sakai et al., 2006).

To obtain gene expression profiles of the study species over a flowering period, we collected leaves and buds sampled from the top canopies of the two S. curtisii trees around noon (Table S1). The concurrent use of both tissues is essential to understand the gene expression dynamics between them, as both tissues are involved in the production and transportation of floral regulatory proteins. These samples, which are part of the collection described in Yeoh et al. (2017), consisted of eight time points (TPs): 08 July 2013 (TP-A), 17 October 2013 (TP-B), 18 December 2013 (TP-C), 14 February 2014 (TP-D), 21 March 2014 (TP-E), 02 April 2014 (TP-F), 14 May 2014 (TP-G), and 17 June 2014 (TP-H). The TPs were selected to represent different flowering stages of the trees, which are categorized based on morphological characters of the buds observed during sampling. In this study, the trees were classified into four stages: (i) vegetative stage, characterized by vegetative buds that would differentiate into leaves, (ii) inflorescence stage, indicated by the presence of inflorescence, (iii) reproductive/flowering stage, marked by the presence of some reproductive buds that would differentiate into flowers, and (iv) abortion/fruiting stage, in which the flowers either developed into fruits or fell off (Fig. S3 and Table S1). Multiple time points were selected for the vegetative stage to capture changes in gene expression leading up to the flowering period. After collection, the samples were immediately submerged in RNAlater solution (Ambion, Austin, TX, USA) and stored at −80 °C until RNA extraction. Using RNASeqPower package (Hart et al., 2013) in R v4.0.3 (R Core Team, 2018), the statistical power of the experimental design was calculated to be 0.38 for a sample size of two (the number of biological replicates used in the study).

Transcriptome sequencing and annotation

Total RNA was extracted from leaf and bud samples using CTAB method described by Kobayashi et al. (2013). The extracted RNA was used to create paired-end cDNA libraries as per manufacturer’s protocol (Illumina, Foster City, CA, USA) followed by RNA-seq using Illumina HiSeq 4000 Sequencer (Illumina). Low quality bases and traces of Illumina adapters were removed using Trimmomatic v0.36 (Bolger, Lohse & Usadel, 2014) prior to de novo assembly using Trinity v2.8.5 (Grabherr et al., 2011). To minimize transcript redundancy, only transcripts with complete open reading frames (ORFs) as predicted by TransDecoder v5.5.0 (http://transdecoder.github.io) were retained. The number of transcripts was further reduced by clustering highly similar protein sequences using CD-HIT v4.8.1 (Fu et al., 2012) (parameters: -n 5 -c 1). The completeness of the assembly was assessed using BUSCO v5.2.1 (Manni et al., 2021). Unless otherwise stated, all bioinformatic analyses in this study were conducted using default parameters.

The non-redundant transcripts, henceforth referred to as unigenes, were annotated by querying against the proteome of A. thaliana (Cheng et al., 2017) using BLASTx v2.10 (Camacho et al., 2009) with E-value cut-off: 1E–10. The annotated unigenes were compared to the A. thaliana flowering genes database (Bouché et al., 2016) to identify their homologs in S. curtisii. Translated unigenes were searched against UniProt database using BLASTp v2.10 (Camacho et al., 2009) with E-value cut-off: 1E–05, as well as Pfam v32.0 (El-Gebali et al., 2018) and PANTHER v14.1 (Thomas et al., 2003) databases using InterProScan v5.36 (Jones et al., 2014). The unigenes were assigned Gene Ontology (GO) terms (Ashburner et al., 2000) associated with the earlier BLASTx matches and visualized using WEGO online tool (Ye et al., 2018). The unigenes were also queried against KEGG database (Kanehisa & Goto, 2000) using BlastKOALA online tool (Kanehisa, Sato & Morishima, 2016).

Differential expression analysis and enrichment tests

Transcript quantification was conducted using Salmon v1.5.1 (Patro et al., 2016). To summarize the expression data, we performed a principal component analysis (PCA) using pcaExplorer package v2.22.0 (Marini & Binder, 2019) in R, which normalized the expression values using DESeq method (Anders & Huber, 2010) prior to PCA. Subsequent bioinformatic analyses were performed on leaf and bud samples separately. We conducted differential expression analysis using DESeq2 software v1.32.0 (Love, Huber & Anders, 2014) by comparing all possible pairwise TPs (Table S2) to identify DEUs (defined as transcripts with absolute log₂ fold change ≥ 1 and false discovery rate (FDR) <1E–03). To identify significantly enriched GO terms (P < 0.05) in the DEUs, enrichment analysis was performed using topGO package v2.40.0 (Alexa & Rahnenfuhrer, 2020) in R. Enriched KEGG pathways (FDR <0.05) were identified using an online tool, KO-Based Annotation System (KOBAS) v3.0 (Bu et al., 2021). To validate the expression profiles of the unigenes, the expression of key flowering gene, FT, homologs was compared to relative expression values measured by qRT-PCR in a previous study in S. curtisii (Yeoh et al., 2017). We utilized ggpubr package in R to perform Spearman’s correlation test and to plot linear regression scatter plots.

Clustering and characterization of differentially expressed unigenes

To summarize the expression profiles of the DEUs in both leaf and bud tissues, the unigenes were clustered based on their expression patterns using coseq package (Rau & Maugis-Rabusseau, 2018) in R. The expression values were normalized using Trimmed Mean of M-values (TMM) method (Robinson & Oshlack, 2010), transformed using centered log ratio (Aitchison, 1982), and clustered using K-mean algorithm (MacQueen, 1967).

For the characterization of the DEUs, we followed the method used in the transcriptome study of S. beccariana (Kobayashi et al., 2013), which compared DEUs to A. thaliana transcriptome data and assumed that the regulation of homologs in both species is similar. This approach enabled functional annotations of DEUs with categories that are not covered by the GO database, such as “response to prolonged moderate drought” (Kobayashi et al., 2013). We compared the DEUs in S. curtisii to 35 gene sets of up- and downregulated genes in A. thaliana that were exposed to various external and endogenous factors (Cao et al., 2006; Gould et al., 2006; Harb et al., 2010; Ma & Bohnert, 2007; Misson et al., 2005; Nemhauser, Hong & Chory, 2006; Osuna et al., 2007; Scheible et al., 2004; see Table S3 for details). Fisher’s exact tests with Bonferroni multiple testing correction were performed to determine the significance in the number of overlapping genes between DEUs and the A. thaliana gene sets (Fig. S4). When the number of DEUs that overlaps with the genes in the gene set is significant, we assumed that the S. curtisii trees used in this study experienced similar conditions as the A. thaliana from which the gene set was obtained.

To determine if the significantly tested A. thaliana gene sets showed specific expression patterns, we further compared the gene sets to DEU clusters using gene enrichment tests (Fig. S5 ; Tables S16–S23, S24–S30). The enrichment tests were performed using Fisher’s exact tests with Bonferroni corrections. We assumed enriched gene sets had the same expression pattern as their significantly associated DEU clusters. Finally, the DEU clusters with significantly enriched gene sets were compared to the temperature and rainfall records to determine potential association between the transcriptome and meteorological data.

RNA sequencing and transcriptome assembly

The transcriptome profile of S. curtisii at different stages of flowering (vegetative, inflorescence, flowering, and fruiting/abortion) was captured by sequencing total RNA from leaves and buds at eight time points (TPs). The sequencing generated 1,595,323,292 reads with an average length of 150 bp which were assembled into 664,053 transcripts. After removing redundant transcripts, the resulting transcriptome assembly comprises 64,219 unigenes with an average length of 2,148 bp and N50 of 2,624 bp (Fig. S6). The GC content of the assembled transcriptome was 43.02%. The BUSCO analysis indicated an almost complete transcriptome assembly in which 95.5% of highly conserved proteins in the Embryophyta lineage were recovered (Table S4). Of the total unigenes, 58,722 (91.44%) were annotated with at least one of the public protein databases (Tabes S5–S7; Fig. S7).

We conducted PCA on normalized unigene expression data to obtain an overview of the variance in gene expression among samples. The analysis clustered the samples based on tissue type, individual tree, and TPs (Fig. 1A). Most of the variance (PC1: 34.71%) can be attributed to transcriptional differences between buds and leaves (Fig. 1A). Variation in transcriptional profiles between the two biological replicates was also high (PC2: 18.63%; Fig. 1A). In bud, 11.29% of the variation (PC3) distinguished majority of the bud samples at vegetative stages (TP-A–D) from the remaining bud samples (TP-E–H; Fig. 1B). The clustering of TP-E with TPs at the inflorescence (TP-F) and later stages (TP-G–H) suggests that the switch to inflorescence may have already begun although morphologically the buds at TP-E are still in vegetative stage. Interestingly, a number of the top contributing unigenes for PC3 (Table S10) are homologs of stress-responsive genes such as EXPANSIN A15 (Wieczorek et al., 2006) and WRKY33 (Zheng et al., 2006).

Differential gene expression and clusters of differentially expressed unigenes

To identify candidate unigenes that are involved in floral regulation of S. curtisii, we conducted a differential expression analysis for pairwise comparison of all TPs. A total of 1,730 and 3,559 DEUs were identified in leaf and bud samples, respectively, with 789 unigenes shared between the two tissues (Tables S11 and S12). The differences in the number of DEUs between the tissues may be due to the regulation of organ-specific unigenes that function in floral morphogenesis, as previously reported in A. thaliana (Smaczniak et al., 2017). To summarize the expression profiles of the DEUs, K-mean clustering was performed, resulting in the identification of six DEU clusters with distinctive expression profiles in leaf transcriptomes (Fig. 2) and five DEU clusters in bud transcriptomes (Fig. 3).

Differential expression of putative flowering-time unigenes

To gain insight into the floral regulatory pathways of S. curtisii, we first identified homologs of A. thaliana flowering-time genes from the annotated unigenes. Of the 51,654 unigenes that could be mapped to A. thaliana genes, 918 unigenes were homologous to 228 non-redundant flowering-time genes (Table 1 and Table S13). This indicates that >75% (228/306) of A. thaliana flowering-time genes have homologs in S. curtisii, similar to the findings in a closely related species, S. beccariana (Kobayashi et al., 2013). Out of the total A. thaliana flowering-time homologs, 103 were differentially expressed throughout the study period (Table S13).

There are 274 flowering-time homologs associated with the photoperiod and circadian clock pathway (Table 1 and Table S13). Photoperiod-associated homologs accounted for over half (15/26 DEUs) of flowering-time DEUs in leaf, while only one third (37/96 DEUs) in bud were from this pathway (Table S13). The number of homologs related to the circadian clock and photoperiod pathway in this study is consistent with our earlier study in S. curtisii leaves despite using fewer time points in the previous study (Suhaimi et al., 2023). We also identified 99 homologs belonging to the vernalization pathway, 34 homologs from the autonomous pathway, and 16 ambient temperature-responsive homologs in the transcriptome. Of the homologs from the vernalization pathway, eight were differentially expressed only in bud, while two homologs of VERNALIZATION 5 (VRN5) were differentially regulated in both leaf and bud (Table 1). Only two homologs from the autonomous pathway were differentially expressed in bud (Table 1). Of the homologs associated with ambient temperature, two homologs of RELATED TO ABI3/VP1 2 (RAV2) were differentially expressed in both leaf and bud (Table S13). Additionally, we identified 41 homologs of integrators in this study, ten of which were differentially expressed only in bud. Two homologs of FT, namely ScFT1 and ScFT2, were differentially regulated in both tissues (Table S13). The expression of these FT homologs is congruent with qRT-PCR data from a previous study (Yeoh et al., 2017; Fig. S7), thus, validating the RNA-Seq results.

Characterization of the differentially expressed unigenes through A. thaliana gene sets enrichment analyses

Only 2.7% of the total DEUs (122/4,500) are homologs of A. thaliana flowering-time genes (Tables S11–13), which suggests that most of them have functions that are not directly related to known flowering mechanisms in the model species. To characterize these DEUs, enrichment analyses were performed on 35 gene sets of A. thaliana (Table 2) and the DEUs separately for leaf and bud samples. We found 25 A. thaliana gene sets to be significantly enriched, that is, these gene sets share significant number of overlapping genes with the DEUs for both S. curtisii leaf and bud (Table 2; Tables S14 and S15). These gene sets include responses to drought and temperature, as well as changes in level of hormones such as jasmonic acid (JA) and gibberellin (GA). We assumed that the trees experienced similar conditions as the A. thaliana from which the gene sets were obtained when the tests were significant.

Expression profiles of the enriched gene sets

To analyze the expression changes of the enriched gene sets across the study period, enrichment analysis was conducted on the gene sets that were found to be significant (Table 2) and each DEU cluster in leaf and bud of S. curtisii (Figs. 2 and 3). We assumed that the expression pattern of gene sets enriched in a particular DEU cluster is represented by that of the cluster. Overall, 12 gene sets were enriched in at least one DEU cluster in leaf (Fig. 2, Table S23) while 13 gene sets were overrepresented in at least one DEU cluster in bud (Fig. 3 and Table S30). Notably, a number of gene sets related to nutrient levels, biotic and abiotic stresses, and hormonal treatments were enriched in DEU clusters that showed similar expression in both trees (Figs. 2 and 3).

Our analysis of DEU clusters in leaf has found that the homologs of A. thaliana genes that are upregulated under medium-term phosphate-limited condition were enriched in DEU cluster 1 (Fig. 2A). The expression profile of this cluster indicates a decrease in phosphorus level between TP-B and TP-D, followed by another decrease prior to flowering. DEU cluster 3 in leaf was overrepresented with homologs of genes commonly upregulated in response to biotic and abiotic stresses, as well as genes upregulated with hormonal treatment, specifically abscisic acid (ABA), indole-3-acetic acid (IAA, or commonly known as auxin), and JA (Fig. 2C). These phytohormones are also involved in plant stress response (Yang et al., 2019). Furthermore, GO terms related to stress response such as “response to reactive oxygen species” and “MAPK cascade”, were enriched in this cluster (Table S34). The expression of the DEUs was elevated from TP-D to TP-F (Fig. 2C), suggesting that the plants were under stress during this period. DEU cluster 4 in leaf was enriched in genes that are upregulated under carbon-limited condition and downregulated with addition of sucrose (Fig. 2D). The expression profile of this cluster suggests that the concentration of carbon in the leaf was decreasing at the onset of flowering, whilst sucrose level was increasing prior to floral induction, specifically at TP-E in C1 and TP-D in C2. Additionally, the KEGG pathway analysis identified the “carbon metabolism” pathway as being enriched in this cluster (Table S34). Finally, DEU cluster 6 in leaf was overrepresented with gene sets responding to drought and changes in the levels of GA and carbon (Fig. 2F). However, the expression patterns of the DEUs in the cluster varied between the two individuals, implying that the response observed was specific to each individual. As a result, the enriched gene sets in the cluster were not further analyzed.

Examination of DEU clusters in bud identified that the homologs of A. thaliana genes upregulated in response to carbon-limited condition were enriched in DEU cluster 1 (Fig. 3A). The expression pattern of the cluster suggests that the level of carbon in bud was depleted as flowering approached (Fig. 3A), which is consistent with the results of enrichment tests in leaf (Fig. 2D). GO terms related to catabolism such as “polysaccharide catabolic process”, “cellular amide catabolic process”, and “allantoin catabolic process”, were also found to be overrepresented in the cluster (Table S35). This aligns with previous research in A. thaliana which demonstrated that carbon starvation leads to an upregulation of genes involved in catabolism (Osuna et al., 2007). Among the DEU clusters in bud, clusters 2 and 4 showed inverse expression patterns in both individuals (Figs. 3B and 3D). Cluster 2 was enriched with homologs of genes that are downregulated under moderate drought, while cluster 4 was not. The expression profile of DEU cluster 2 suggests that the Shorea individuals were experiencing mild drought between TP-D and TP-F. Additionally, both clusters 2 and 4 were enriched with homologs of genes that respond to elevated level of GA. The expression profiles of these two clusters indicate an increase in GA level between TP-D and TP-F. Both clusters were also enriched with GO terms related to stress, signaling, growth, and floral regulation, such as “response to reactive oxygen species”, “vasculature development”, “MAPK cascade”, and “signal transduction” in cluster 2, as well as “response to hypoxia”, “response to light stimulus”, “microtubule cytoskeleton organization”, and “cytokinesis” in cluster 4 (Table S35). Our PCA also revealed that the majority of top contributing unigenes to PC3 belong to these two clusters (Table S10), highlighting the importance of the unigenes in these clusters in the transition from vegetative to reproductive stages. In DEU cluster 3, we found homologs of A. thaliana genes that are typically upregulated in response to diverse abiotic and biotic stresses, as well as genes that are affected by cytokinin and JA treatment (Fig. 3C). However, this cluster showed inconsistent expression profiles between the two S. curtisii individuals (Fig. 3C), hence, the enriched gene sets in this cluster were not examined further. Lastly, DEU cluster 5 in bud was overrepresented with homologs of genes upregulated under moderate drought and downregulated with increasing temperature (Fig. 3E). The expression profile of this cluster indicates an increase in temperature between TP-D and TP-E (Fig. 3E), which is in line with the enrichment of relevant GO terms such as “response to stress” and “response to stimulus” (Table S35).

Association between meteorological data and transcriptome profiles

Prolonged moderate drought and drop in temperature have been hypothesized as proximate cues for GF (Ashton, Givnish & Appanah, 1988; Brearley et al., 2007; Chen et al., 2018; Yeoh et al., 2017). Therefore, we examined meteorological data to look for these signals prior to flowering. These data were then compared to the expression patterns of DEU clusters that were significantly enriched with A. thaliana gene sets related to temperature or drought (Figs. 2 and 3) to infer the regulatory pathways that respond to these climatic cues.

The daily mean temperature at the nearest meteorological station situated 11 km from the study site was 26.8 ° C, while the total annual precipitation at the nearest hydrological station 6 km from the study site was 2,272 mm. The lowest minimum daily temperature (18.3 ° C) was recorded in early February around TP-D (Fig. 4B). We also calculated 30-day cumulative rainfall to identify drought periods, which are defined as 30-day cumulative rainfall of below 40 mm (Sakai et al., 2006). The cumulative rainfall data showed a period of drought beginning before TP-D and lasting until TP-E (Fig. 4C). Therefore, the expression patterns of leaf DEU cluster 3 (Fig. 2C) as well as bud DEU clusters 2, 4, and 5 (Figs. 3B, 3D and 3E) were examined.

Among these DEU clusters, only DEU cluster 5 in bud was enriched with homologs of genes downregulated with increased temperature (Fig. 3E) and its expression pattern implies that the plants were responding to rising temperature around TP-D. This pattern concurs with the daily mean temperature data, which showed a slight increase within the same period (Fig. 4A). Although the increase in temperature was subtle, it was not unexpected for the plants to be able to capture the slight change in temperature stimuli as it has been suggested that trees in tropical regions are more sensitive to change in their surroundings (Singh & Kushwaha, 2016). Of the DEU clusters overrepresented with homologs of genes responding to drought, only DEU cluster 3 in leaf (Fig. 2C) and cluster 2 in bud (Fig. 3B) displayed expression patterns consistent with the drought period (Fig. 4C). This indicates that the trees were experiencing a moderate drought between TP-D and TP-E. The expression pattern of DEU cluster 5 in bud (Fig. 3E), on the other hand, did not match the cumulative rainfall data.

In essence, we assembled a transcriptome with majority of the unigenes annotated. The examination of these annotated unigenes, along with the differential expression analysis, identified homologs of flowering-time genes associated with various pathways, such as photoperiod, circadian clock, and ambient temperature responses. Subsequent enrichment analyses and the expression profiles of DEU clusters suggested that the plants responded to factors such as phosphate and carbon availability, and hormonal dynamics. Additionally, the examination of the meteorological data in relation to the DEU clusters indicated that the trees exhibited changes in transcriptional expression in response to rising temperatures and a brief spell of drought.

Proximate environmental cues for general flowering

The increasing rate of climate change, characterized by higher temperature fluctuations (IPCC, 2022) and more frequent and severe droughts (Dai, 2013; Trenberth et al., 2014), among others, is expected to have a disproportionate impact on the phenology of plants and animals in tropical regions (Santer et al., 2018). The fitness and abundance of tropical communities, including GF dipterocarp species, can be adversely affected by a slight temperature increase (Lister & Garcia, 2018; Numata et al., 2022). As these plant species is a valuable source of sustenance for a variety of seed and fruit predators including wild pigs and rodents (Ickes, 2001; Miura, Yasuda & Ratnam, 1997), alteration in their phenology could have potential cascading effects on the tropical forest communities.

Long term phenological records of tropical rainforests in the Malay Peninsula have shown that flowering events are always preceded by a drop in daily minimum temperature after a brief period of drought (Ashton, Givnish & Appanah, 1988; Sakai et al., 2006). In fact, mathematical models that take into account the synergistic effects of both drought and a drop in daily minimum temperature have been found to be reliable at predicting GF occurrences (Chen et al., 2018; Yeoh et al., 2017). The cumulative rainfall data (Fig. 4C), along with the expression profile of DEU clusters that were enriched with homologs of drought-responsive genes (Figs. 2C and 3B) indicate that the S. curtisii in our study experienced water stress due to a moderate drought lasting from before TP-D until TP-E. The upregulation of ScFT1 (Unigene: TRINITY_DN6921_c0_g2_i2) in leaf shortly after drought began (TP-D; Table S11) and the upregulation of ScFT1 in bud post-drought (TP-E; Table S12) further corroborated the hypothesis that drought is one of the floral triggers in dipterocarps (Kobayashi et al., 2013; Sakai et al., 2006; Yeoh et al., 2017). Moreover, the observed expression changes in homologs associated with elevated levels of ABA, IAA, and JA hormones (Fig. 2C) that coincides with the occurrence of drought also supported the aforementioned hypothesis. Roles of these hormones in the floral regulation and signaling pathways of plants responding to stress, including drought, have been established in many other species such as A. thaliana, rice and Citrus spp. (reviewed in Singh & Laxmi, 2015; Yang et al., 2019). Among these hormones, ABA is particularly important in the hormonal response to drought stress (Gupta, Rico-Medina & Caño Delgado, 2020) as it regulates stomatal opening to reduce water loss (Okamoto et al., 2013). Interestingly, ABA has also been reported to function in drought-induced flowering in A. thaliana (Riboni et al., 2016) and Citrus spp. (Khan et al., 2022; Li et al., 2017). Auxin (IAA) plays a role in modulating leaf water uptake and regulating the expression of antioxidant enzymes to detoxify reactive oxygen species (Shi et al., 2014). Similarly, JA contributes to abiotic stress responses by increasing resistance to oxidative stress (Wu et al., 2012) and maintaining downstream drought-responsive pathways such as those related to modulation of shoot biomass and photosynthetic rate (Wang et al., 2021).

Despite a sharp drop in daily minimum temperature observed prior to TP-D (Fig. 4B), we did not find any gene set associated with differential regulation in response to decrease in temperature to be overrepresented in any of the DEU clusters (Figs. 2 and 3). However, several studies have shown that temperature drops always precede GF (Ashton, Givnish & Appanah, 1988; Numata et al., 2003; Yasuda et al., 1999), and that the occurrence of drought alone does not always result in GF (Chen et al., 2018; Yeoh et al., 2017). Therefore, it is possible that S. curtisii employs different temperature-responsive mechanisms than A. thaliana. This scenario is not unlikely, as similar divergence has been reported in other floral regulatory pathways. For example, FLOWERING LOCUS C (FLC) encodes a key floral repressor in the vernalization pathway of A. thaliana, grasses have developed vernalization response that is completely independent of FLC, relying instead on four central genes: VERNALIZATION 1–3 (VRN1–3) and VEGETATIVE TO REPRODUCTIVE TRANSITION 2 (VRT2) (reviewed in Trevaskis et al., 2007). Further studies using transcriptome profiles of GF dipterocarp trees treated with low temperature and drought under controlled environments could help to verify the potential of low temperature as a flowering cue and further prove that a short period of drought can induce flowering in the species. However, this would be a practically challenging and costly endeavor for tall tropical trees that flower irregularly.

Additionally, several homologs of flowering-time genes were found to be differentially expressed after the short dry spell (Figs. 2 and 3; Tables S11–13). Given the accumulating evidence of organisms and populations adapting to climate change through phenotypic plasticity and genotypic evolution (reviewed in Peñuelas et al., 2013), further research focusing on these homologs could help us understand how these rapid environmental changes will impact the reproductive phenology of GF dipterocarp species and how this phenomenon will evolve in the future.

Resource dynamics in Shorea during flowering time

Accumulated nutrients have been proposed as a limiting factor in GF, which could explain why some mature trees do not flower during GF (Ichie et al., 2013; Kelly & Sork, 2002). In the present study, we found molecular evidence showing changes in the expression levels of homologs associated with a decrease in phosphorus level prior to flowering (Fig. 2A). This finding is consistent with a previous report that investigated mineral nutrient storage dynamics prior to mast reproduction in another dipterocarp species, Dryobalapnos aromatica (Ichie & Nakagawa, 2013). Although we identified homologs of A. thaliana genes responsive to nitrogen to be overrepresented in the DEUs of S. curtisii (Table 2), these unigenes did not show any specific expression pattern (Figs. 2 and 3). In addition to phosphorus and nitrogen, carbohydrates have been demonstrated to be utilized for floral initiation and subsequent development of floral organs (Peng & Iwahori, 1994). The results of our cluster analysis indicate changes in the expression level of unigenes correlated with an increase in sucrose level in leaves during the brief period of drought (Fig. 2D). An elevated level of sucrose in leaves during drought have been reported in other trees such as East Asian white birch (Bhusal et al., 2021) and tropical starfruit (Pingping, Chubin & Biyan, 2017). Furthermore, the cluster analysis identified changes in expression levels of homologs correlated with a reduction in carbon concentrations in bud prior to flowering time (Fig. 3A). Previous studies have suggested that flower bud act as carbohydrate ‘utilizing sinks’, metabolizing rather than storing carbohydrate (Eshghi et al., 2007; Ho, 1988).

In the current study, we inferred the fluctuation of nutrient levels based on associated or correlated transcript expressions. To gain a more comprehensive understanding of the roles and dynamics of plant resources during flowering time, future studies should consider integrating nutritional content information such as nitrogen, phosphorus and carbon concentrations. The collection of sufficient tissues from carbohydrate sinks such as stems and trunks (Kozlowski, 1992) without damaging the branches was not feasible in this study due to the need for frequent repeated sampling to capture the floral initiation time points of these trees with unpredictable flowering intervals. Furthermore, canopy access for tall tropical trees has thus far been limited, hence observations and sample collections in our study were done by climbing approximately 40 m above the ground canopy. It is also important to acknowledge the low statistical power of this study due to the limited number of biological replicates and the absence of technical replicates. In addition, due to logistical constraints, rainfall measurements in the current study were obtained from the nearest hydrological station located adjacent to the forest reserve, rather than being recorded at the specific location of each individual tree. Therefore, to address these limitations, it would be essential to employ improved monitoring and sampling techniques e.g., by utilizing technology such as unmanned aerial vehicles (drones) for sample collection, to facilitate future research. Notwithstanding these challenges and limitations, we were able to conduct the present study on two individual trees and obtain clear and consistent results.

Based on the findings in the current study and earlier reports (Chen et al., 2018; Ichie et al., 2013; Kobayashi et al., 2013; Yang et al., 2019; Yeoh et al., 2017), we propose a preliminary framework that depicts the interactions between external environmental cues, phytohormones, nutrient resources and expression of florigen, FT in Shorea (Fig. S9) during GF. According to this framework, flowering signals such as drought (and potentially, drops in daily minimum temperature) accumulate over a period of 2–3 months prior to floral initiation. During this period, the levels of stress-associated phytohormones, including IAA, ABA, and JA, increase, leading to the upregulation of the key flowering gene, FT in the leaf. Subsequently, FT expression in the bud is elevated, thus initiating the transition to the reproductive stage. Further research is needed to develop and complete this proposed framework and determine the extent of its applicability to other GF dipterocarp species.