Characterization of leaf transcriptome in a tropical tree species, Shorea curtisii, over a flowering season

General flowering (GF) is a synchronous flowering event in the Southeast Asian tropical rainforests that occurs at irregular intervals of multiple years. The unpredictable intervals of GF raise conservation concerns for these under-researched forests with rich economically and ecologically important species. In this study, the leaf transcriptome of a GF species, Shorea curtisii obtained from three time points – before and after floral initiation, and post flowering stage – was sequenced. We assembled 243,759,478 sequencing reads into 39,943 non-redundant unigenes including 677 putative homologs of Arabidopsis thaliana flowering-related genes. Differential expression analysis conducted on pairwise comparisons of the time points identified 930 differentially expressed unigenes, which includes 17 flowering-related homologs. The differential expression of unigenes with significant enrichments of functions related to drought corroborated the involvement of drought as an environmental cue for GF. The outcomes of this study offer an insight into the conservation of floral regulatory genes and pathways in Shorea and could be used as a model to better understand the floral initiation cues and regulation of GF trees.

Methods
Leaf samples at top layer of canopy from two S. curtisii individuals, C1 and C2 were collected from Semangkok Forest Reserve (2°58’N, 102°18’E). We selected three time points (TPs) corresponding to before and after floral initiation, and abortion/fruiting stage. Samples were collected at each TP around midday, soaked in RNAlater reagent (Ambion, USA), and stored at –80°C.

Total RNA was extracted from leaf samples using the method described by Kobayashi et al. (2013). RNA-seq was performed using Illumina HiSeq 4000 Sequencer (Illumina, USA). After removal of low quality bases by Trimmomatic v0.36 (Bolger et al., 2014), the sequence reads from each sample were combined and de novo assembled using Trinity v2.8.5 (Grabherr et al., 2011). Only non-redundant transcripts with complete open reading frame were retained by using TransDecoder v5.5.0 (http://transdecoder.github.io) and CD-HIT v4.8.1 (Fu et al., 2012). The quality of the assembly was assessed with BUSCO v5.2.1 (Simao et al., 2015).

The unigenes were queried 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 searched against A. thaliana flowering genes database (Bouché et al., 2016) to identify their homologs in S. curtisii. Translated unigenes were queried against other public protein databases using BLASTp v2.10 (Camacho et al., 2009) and InterProScan v5.36 (Jones et al., 2014). We also searched the unigenes against Gene Ontology (GO; Ashburner et al., 2000) and KEGG pathways (Kanehisa & Goto, 2000) to characterize the unigenes. GO classification is divided into cellular component (CC), molecular function (MF), and biological process (BP) categories.

The unigenes were quantified by Salmon v1.5.1 (Patro et al., 2016) and subjected to differential expression analysis using DESeq2 v1.32.0 (Love et al., 2014). The samples were compared in a pairwise manner and differentially expressed unigenes (DEUs) with absolute log2 fold change ≥ 1 and FDR < 0.05 were identified. Significantly enriched GO terms (P < 0.05) and KEGG pathways (FDR < 0.05) in the DEUs were also identified using topGO package v2.40.0 (Alexa & Rahnenfuhrer, 2020) and KOBAS v3.0 (Bu et al., 2021), respectively.