Genome-wide transcriptome analysis reveals the molecular mechanism of high temperature-induced floral abortion in Litchi chinensis
BMC Genomics volume 20, Article number: 127 (2019)
Warm winter and hot spring attributed to global warming affected floral development and may induce floral abortion, resulted in poor flowering in litchi. To identify genes potentially involved in litchi floral abortion, six RNA-sequencing (RNA-Seq) libraries of the developing panicles (DPs) under low temperature (LT) conditions and the shrinking panicles (SPs) under high temperature (HT) conditions were constructed.
3.07–8.97 × 106 clean reads were generated. Digital expression of the DPs with that of the SPs was compared. As a result, 1320 up-regulated and 981 down-regulated differentially expressed genes (DEGs) were identified. From the enriched GO-term, 54 temperature responsive DEGs, 23 hormone homeostasis- or biosynthesis-related DEGs, 137 hormone signal transduction or responsive DEGs, and 18 flowering-related DEGs were identified. Partial Least Squares Structural Equation Modeling (PLS-SEM) analysis indicated that the effects of hormone-related DEGs on NACs, MYBs, WRKYs were stronger than that on flowering-related DEGs. Expression pattern analysis of the inflorescence in ‘Nuomici’ and ‘Huaizhi’ under LT and HT conditions showed that genes homologous to AIL6 (LcAIL6), LHY (LcLHY), MED16 (LcMED16), SKIP20 (LcSKIP20), POD20 (LcPOD20) in the two cultivars had similar expression trends.
This study elucidated the transcriptome in the HT-induced floral abortion and identified key genes involved in the process. NACs, MYBs, WRKYs may act as central players involved in the HT-induced floral abortion underlying hormonal control. Increased transcript in LcLHY, LcMED16, LcSKIP20, LcPOD20 and decreased transcript in LcAIL6 might be related to the inhibition of floral development. Our studies provided potential genes for the future molecular breeding of new cultivars that can reduce floral abortion under warm climates, and a novel clue to reveal the relationship of biological processes based on the RNA-seq data using PLS-SEM.
Litchi is an evergreen woody tree cultivated in subtropical and tropical regions and produces arillate fruits with sweet, translucent, juicy flesh . However, unstable flowering causing irregular bearing affects the litchi industry. Up till now, many studies indicated that chilling is an irreplaceable factor for litchi flowering [2, 3]. In winter, litchi trees are subjected to low temperature for floral induction. After winter temperature treatment, litchi floral buds which are a mix of axillary and apical panicle primordia, leaf primordia or rudimentary leaves, break and the panicle primordia start to develop . At this time point, whether panicle primordia can develop continually depends on environmental temperature. If the temperature is still low enough, the panicle primordia develop, and the rudimentary leaves cease growing and abscise automatically. The mixed buds will develop into pure/high-quality panicles . However, warm winter and hot spring attributed to global warming frequently happen. The rudimentary leaves may develop into fully expanded leaves (panicle leaves), and the panicle primordia stop developing and shrink . These mixed buds may become leafy panicles or even vegetative shoots, resulting in poor flowering . How to control panicle leaf growth and inhibit floral abortion is of great importance to litchi production. In practice, ethephon is normally used to control panicle leaf growth under high-temperature conditions . Methyl viologen dichloride hydrate, a reactive oxygen species (ROS) producer, and sodium nitroprusside, a nitric oxide (NO) donor are also proved to effectively inhibit panicle leaf growth and floral abortion . Therefore, understanding the regulatory mechanism of floral abortion and panicle leaf controlled by low temperature or by chemicals is of great importance for flowering regulation in litchi.
We have previously focused on the feature of the low temperature-, ROS-, and NO-induced senescence of panicle leaves, and found that programmed cell death (PCD) is involved in the panicle leaf growth regulation by ROS and NO, as well as by low temperature in litchi [6, 8]. We also screened PCD related genes from our RNA-seq data sets of the ROS-treated rudimentary leaves , and identified a litchi homolog MCII (LcMCII-1). Silencing LcMCII-1 by virus-induced gene silencing (VIGS) delayed ROS-dependent senescence. Ectopic over-expression of LcMCII-1 in transgenic Arabidopsis promoted ROS-dependent and natural senescence , suggesting that LcMCII-1 is positively involved in the regulation of rudimentary leaf senescence in litchi. These results provide a new target for the future molecular breeding of new cultivars that can produce high-quality flowers in warmer climates. However, other than studies on panicle leaf control, how the panicle primordia abort or stop developing under high-temperature conditions is still poorly known.
RNA sequencing (RNA-seq) technology is a powerful tool for plant studies. It has been used to identify key candidate genes in response to vernalization of oriental lily, and MADS-box family genes related to organ development and stress resistance in Brassica rapa [11, 12], compare transcriptome between low- and high-cadmium-accumulating Brassica chinensis genotypes in response to cadmium stress , reveal co-expression networks of transcription factors and wood component genes in Populus trichocarpa .
In this study, we focused on the molecular mechanism of floral abortion under high-temperature conditions. We used RNA-Seq technology to compare transcriptomes between the shrinking panicles (SPs) and the developing panicles (DPs). For RNA-Seq library construction, we used ‘Nuomici’ litchi as it is sensitive to HT in which floral buds easily cease to develop, and the panicle leaves growth vigorously according to our previous studies [4, 6], so that we could get more aborted inflorescences for RNA-seq. The litchi potted trees of ‘Nuomici’ were transferred to growth chambers for low temperature (LT) treatment to promote panicle/floral development and high-temperature treatment (HT) to induce panicle/floral abortion. The expression of the DPs and SPs was compared to identify differentially expressed genes (DEGs) potentially involved in litchi floral abortion. Then the expression pattern of the candidate genes in different cultivars under HT or LT conditions were determined. Partial Least Squares Structural Equation Modeling (PLS-SEM) analysis was also performed to quantify the relationship among the biological processes based on the transcript levels of the DEGs. We aim to reveal the molecular mechanism of the floral abortion, and to provide potential genes for the future molecular breeding of new cultivars that can set fruit in warmer climates.
Morphology of DPs and SPs
To identify genes involved in floral abortion, ‘Nuomici’ trees were grown under HT to induce floral abortion and to produce SPs, or LT to promote floral development and produce DPs as a control. As a result, the percentage of flowering terminal shoots of the HT-treated trees (50.1%) was significantly lower than that of the LT-treated trees (85.3%) at 0.05 probability level. Figure 1 shows the morphology of a DP in LT and an SP in HT. The DP shows a main-order axis and two second-order axes in the inflorescence. Differentiation of axillary panicles progresses basipetally, showing the activity of the apical meristems, while the SP shows a delay or cease of development, indicating a low activity of meristem.
Digital transcriptome analysis and DEGs identification
Six RNA-seq. libraries of DP and SP were constructed to identify genes potentially involved in HT-induced floral abortion. As shown in Table 1, we have generated 3.12–9.27 × 106 raw reads and got 3.07–8.97 × 106 clean reads from the libraries, with clean read ratios more than 95%. All the clean reads were mapped to the litchi transcriptome with mapping rates more than 55%. As a result, we obtain 1.91–5.24 × 106 mapped reads for the six libraries.
The unique match reads from the 6 libraries were remapped to the reference sequences and the unigene expression was normalized to FPKM. A linear regression analysis result of the fold-change in the gene expression ratios between RNA-seq and qRT-PCR was performed. A significantly positive correlation was found between RNA-seq data and qRT-PCR data, suggesting a reliable transcriptome analysis using RNA-seq (Additional file 1: Figure S1).
Correlation analysis of the DP and SP data suggested high repeatability of the sequencing samples (Additional file 2: Figure S2). DEGs between DP and SP were detected by DESeq2. As shown in Fig. 2, from the DEGs of the DP to those of the SP, 1320 up-regulated genes and 981 down-regulated genes were identified.
GO-term analysis of DEGs
GO-term analysis was performed among the DEGs. The top 10 enriched GO terms are shown in (Additional file 3: Table S1). All the DEGs were classified into the cellular component, biological process, and molecular function. Under the biological process category, ribosome assembly, mature ribosome assembly, and ribonucleo protein complex assembly were the top abundant subcategories. Under cellular component, cytosolic part, cytosolic ribosome, and ribosomal subunit were the top abundant subcategories. For molecular function category, structural constituent of ribosome, structural molecule activity, and copper ion binding were the top abundant subcategories. Ribosome assembly related terms were among the most significantly enriched ones (Additional file 3: Table S1). Intriguingly, as shown in Table 2, several temperature-, hormone-, and flowering-related terms such as response to temperature stimulus, response to hormone, and regulation of floral meristem growth were significantly enriched in the biological process. Given their direct link to the variable we used and the physiological distinctions we observed in our experiment, we focused on those terms in our downstream analyses.
Identification of flowering related genes involved in HT-induced floral abortion
From the enriched flowering related GO-term including regulation of floral meristem growth, regulation of inflorescence meristem growth, floral meristem growth, long-day photoperiodism (Table 2), we identified 8 flowering-related DEGs (Fig. 3). They encode homologous proteins including one Zinc-finger homeodomain protein (ZHD4), two Probable Glutamate Carboxypeptidase 2 (AMP1), one Nuclear Cap-binding Protein Subunit (ABH1), one Mediator of RNA Polymerase II Transcription Subunit 16 (MED16), two Late Elongated Hypocotyls (LHY), and one Cryptochrome-2 (CRY2). We also identified 10 other flowering related genes from the DEGs, such as Agamous-like MADS-box protein AGL8 encoding gene (AGL8), Protein EARLY FLOWERING 4 encoding gene (ELF4), Zinc finger protein CONSTANS-LIKE 12 encoding gene (COL12), and AP2-like ethylene-responsive transcription factor AINTEGUMENTA encoding gene (ANT). Compared to the DP, most of the flowering-related genes in the SP show down-regulated trends (13/18), such as the AINTEGUMENTA-like AIL5 and AIL6, AGL8, COL12, and ANT (Fig. 3).
Identification of the plant hormone-related genes involved in HT-induced floral abortion
From the enriched GO-term related to the control of plant hormone levels, we identified 23 DEGs that might control the hormone levels under different temperature conditions (Fig. 4). These enriched GO-terms are the auxin homeostasis, brassinosteroid metabolic process, steroid biosynthetic process, phytosteroid biosynthetic process, brassinosteroid biosynthetic process, brassinosteroid homeostasis, regulation of salicylic acid biosynthetic process, salicylic acid biosynthetic process (Table 2). The DEGs might be involved in hormone homeostasis, such as the Indole-3-acetic acid-amido synthetase GH3.6 encoding genes (GH3.6) that can prevent free IAA accumulation. They also might be involved in hormone biosynthesis process, such as the Delta (24)-sterol reductase encoding gene DIM, or they might control the hormone biosynthesis process, such as the Lipase-like PAD4 (PAD4) encoding gene (Fig. 4).
From the enriched plant hormone-related GO-term, we identified 137 plant hormone signaling or hormone responsive DEGs (Additional file 4: Figure S3; Additional file 5: Figure S4). These enriched GO-terms are salicylic acid mediated signaling pathway, response to jasmonic acid, response to abscisic acid, response to cytokinin, and response to hormone (Table 2). The DEGs encode proteins such as WRKY transcription factor (WRKY), Transcription factor MYB108 (MYB108), transcription factor MYB32 (MYB32), NAC transcription factor 56 (NAC056), Ethylene-responsive transcription factor (RAP2–3), Homeobox-leucine zipper protein (ATHB-7). Compared to the DP (from DP to SP), most of the hormone-controlled transcription factor encoding DEGs were up-regulated (Additional file 4: Figure S3). On the whole, more up-regulated DEGs were found in the enriched plant hormone-related GO-term (Additional file 4: Figure S3; Additional file 5: Figure S4).
Pathway probably involved in the HT-induced floral abortion underlying hormonal control
To reveal a frame work of the HT-induced floral abortion by the regulation of hormone, we performed PLS-SEM analysis using the transcriptome data of the hormone-related and flowering-related DEGs as described above, and the main hormone-responsive transcription factor encoding DEGs including NACs, MYBs, and WRKYs. We choose PLS for data analysis because it not only works well on removing the collinearity among variables but also works steadily with many kinds of database as well. In PLS-SEM, most of the latent variables were significant, outer model blocks have made an adequate explanation for latent variables (Additional file 6: Table S2). Average variance extracted (AVE, an indicator for converge validity) and composite reliability (indicator for internal consistency reliability) were higher than 0.5 and 0.7 (Additional file 7: Table S3), respectively, indicating that the model was reliable . As shown in Fig. 5, there are two kinds of pathways for each model, the direct pathway and indirect pathway. In indirect pathways, the effects of hormone-related DEGs on the three main transcription factors and the effects of these transcription factors on flowering-related DEGs were positive or negative. All of them had a high factor loading such as 0.945, 0.993, and 0.917 in the CTK model, higher than the second order in indirect pathways (0.130, − 0.359 and − 0.208). Apart from IAA-related DEGs, the CTK-, JA-, SA-, BR-related DEGs were connected with flowering-related DEGs directly. Though IAA-related DEGs were not related to flowering-related DEGs directly, they were related to these DEGs through NACs and MYBs. On the whole, the effects of hormone-related DEGs on NACs, MYBs, WRKYs were stronger than that on flowering-related DEGs.
Identification of temperature responsive genes involved in HT-induced floral abortion
From the enriched temperature-related GO-term including response to temperature stimulus, response to cold, cold acclimation, we identified 54 DEGs potentially involved in the HT-induced floral abortion. They encode homologous proteins such as Delta (8)-fatty-acid Desaturase 2 (SLD2), Cold-regulated 413 Inner Membrane Protein 1 (COR413IM1), Mannose-1-phosphate Guanylyltransferase (CYT1), Heat Stress Transcription Factor (HSFA2), Galactinol Synthase 1 (GOLS1), Dehydration-responsive Element-binding Protein 1A (DREB1A), Peroxidase 3 (PER3). Compared to the DP, most of the temperature-related genes of the SP showed up-regulated trends (33/54) (Fig. 6), such as the PER3, Heat stress transcription factor A-2 (HSFA2), Dehydration-responsive element-binding protein 1A (DREB1A) encoding genes. Some temperature responsive genes showed down-regulated trends, such as the Delta(8)-fatty-acid desaturase 2 (SLD2), Glycine-rich RNA-binding protein (RZ1A), and Probable glutathione S-transferase (HSP26-A).
Expression profiles of the candidate genes in different cultivars under HT and LT conditions
To study detail expression patterns of the candidate genes in DPs and SPs, expression patterns of the candidate genes in the two widely cultivated ‘Nuomici’ and ‘Huaizhi’ under HT and LT conditions were studied. Flowering conditions of ‘Nuomici’ and ‘Huaizhi’ are shown in Table 3. Percentage of flowering terminal shoots in the LT-treated ‘Nuomici’ and ‘Huaizhi’ trees was significantly higher than that of the HT-treated trees. However, under HT condition, the percentage of flowering terminal shoots in ‘Huaizhi’ was higher than that in ‘Nuomici’. As to the percentage of leafy panicles, that of the LT-treated ‘Nuomici’ was the lowest, while that of the HT-treated ‘Nuomici’ was the highest. Compared to LT-treated trees, HT-treated trees of the two cultivars had smaller panicles according to the length of panicles and the length of axillary panicles, indicating that LT inhibited flowering.
Twelve genes identified from the flowering-, hormone-, and temperature response-related DEGs were selected. The expression pattern of the inflorescence in ‘Nuomici’ and ‘Huaizhi’ at the four panicle developmental stages under HT and LT (stages 1–4) according to Yang et al.  were determined. As shown in Figs. 7 and 8, the expression patterns of the genes homologous to AIL6 (LcAIL6) in the two cultivars are similar. Those in LT showed increasing trends, while in HT showed decreasing trends and remained at lower levels, showing positive correlation to flowering. Genes homologous to LHY (LcLHY), MED16 (LcMED16), SKIP20 (LcSKIP20), POD20 (LcPOD20) in the two cultivars showed similar expression trends. Those in HT showed increasing trends, while those in LT showed stable and remained at lower levels, showing a negative correlation to flowering. Similar expression patterns were found in LcCRY2 in ‘Nuomici’, LcWRKY, LcLAP, and LcMYB32 in ‘Huaizhi’. In ‘Nuomici’, floral buds under HT condition had relatively higher levels of LcNAC045 and LcMYB32 expression. However, no matter in ‘Nuomici’ or ‘Huaizhi’, LcNAC100 and LcPOD4 did not show any regular trends. On the whole, 5 out of the 12 candidate genes showed similar expression trends in the two cultivars and were correlated with flowering, and another 5 candidate genes in ‘Nuomici’ or ‘Huaizhi’ showed a similar relationship with flowering.
Warm winter always encourages vegetative growth and repress flowering formation of the evergreen fruit trees in tropical and subtropical areas. Litchi is an evergreen fruit tree whose floral induction and differentiation occur in winter and early spring. There exists competition between reproductive growth and vegetative growth in litchi similar to that of other evergreen fruit trees. Litchi floral buds consist of apical panicle primordia, axillary panicle primordia, and rudimentary leaves. These rudimentary leaves can expand quickly under high-temperature conditions, and then the panicle primordia stop developing and shrink. As a result, the floral buds develop to leafy panicles with few flowers or even reverse to vegetative shoots. Therefore, understanding the regulatory mechanism of floral abortion is as important as that of the transformation from vegetative meristems (VM) to inflorescence meristems (IM) for flowering regulation in litchi. In our previous study, we performed RNA-seq and obtained a global transcriptome of the apical meristem, revealed potential gene networks controlling the transformation from VM into IM in litchi . In the present study, Genome-wide transcriptome analysis of the floral buds of the ‘Nuomici’ litchi trees under low and high-temperature conditions was carried out. We used ‘Nuomici’ as it is sensitive to HT in which floral buds easily cease to develop. Six RNA-Seq libraries were constructed, and 3.07–8.97 × 106 clean reads were generated. To identify genes potentially involved in litchi floral abortion, we compared digital expression of the genes in the DP with that in the SP and identified 1320 up-regulated and 981 down-regulated DEGs. Then GO-term analysis was performed among the DEGs. In the biological process, we found that several temperature-, hormone-, and flowering-related terms were significantly enriched. Mainly basing on these terms, we identified flowering-related, hormone-related, and temperature responsive DEGs that might be involved in floral abortion.
Litchi panicle development undergoes inflorescence and floral organ development. These processes depend on floral meristem activity. The ZHD4, also called FLORAL TRANSITION AT THE MERISTEM2 (FTM2), is expressed in the meristem and correlates with floral meristem growth in Arabidopsis . In this study, microstructure observation of the floral buds indicated that floral meristem growth decreased or ceased under HT condition. Accordingly, the transcript level of the litchi homolog ZHD4 was reduced, suggesting that ZHD4 might be involved in the HT-induced floral abortion. Organ formation in plants is dependent on stem cell niches (SCNs) located in the meristems. Meristems show a functional zonation along the apical-basal axis and the radial axis. Organ primordia are formed in the circular peripheral zone (PZ) from stem cell descendants in which differentiation programs are activated. SCN actively suppresses stem cell identity in the PZ to allow organ development. In Arabidopsis, AMP1 was proven to suppress ectopic SCN formation in the PZ . ANT function in floral initiation and development, mutations in ANT, cause a decrease in floral organ number and alterations in the appearance of all floral organs . ANT acts redundantly with AP2 and function as a class A gene with regard to specification of petal identity . The AINTEGUMENTA-like (AIL) genes are expressed primarily in actively dividing tissues. AIL6 is expressed in inflorescence meristems and flowers . ant ail6 mutants display a delay in the meristem identity transition and LFY induction. ANT and AIL6 transcription factors bind to the LFY promoter to control LFY mRNA accumulation, promote floral meristem identity transition . In our litchi transcriptomic dataset, we identified AMP1, ANT, and AIL6 whose expression was suppressed in HT compared to that in LT. It is suggested that they might be related to the cease of floral development under HT condition.
LHY is an essential clock component that plays an important role in photoperiodic flowering by controlling the rhythmic expression of flowering-time genes in Arabidopsis . It repressed the floral transition under short-day and long-day conditions . Our results showed that expression of the litchi homology LcLHY was higher in HT than that in LT. Accordingly, qRT-PCR analysis showed that LcLHY expression in inflorescence at HT was higher than that at LT in both the ‘Nuomici’ and ‘Huaizhi’. Contrarily, floral differentiation was repressed in HT. It is suggested that the increased expression in LcLHY might be related to the cease of floral development. It seems that floral abortion is regulated by a complicated gene network.
Plant hormones are signal molecules produced in plant cells and are expressed at extremely low concentrations but can regulate the formation of flowers, stems, leaves, control the development and ripening of fruit, and response to biotic and abiotic stresses . Plant hormone levels are related to hormone homeostasis, biosynthesis, and the regulation of the biosynthetic process. In this study, we found that many hormone level regulated genes were differentially expressed between HT and LT conditions. GH3.6 encodes an IAA-amino synthetase that prevents free IAA accumulation, maintains auxin homeostasis by conjugating excess IAA to amino acids . DIM encodes Delta (24)-sterol reductase who plays a critical role in the general process of plant cell elongation, is involved in the conversion of 24-methylenecholesterol to campesterol, an early precursor of brassinolide . PAD4 encodes a lipase-like PAD4 participating in a positive regulatory loop that increases SA levels . They all have higher expression levels in the SPs than in the DPs. In accordance with the gene expression profiles, it is suggested that HT treatment could change hormone homeostasis related gene expression, alter IAA, BR, and SA homeostasis conditions. These hormonal conditions might be related to the cease of floral development.
The plant hormones act as signals are transmitted to the nuclear by series signal transduction components to activate gene expression and result in physiological changes. Transcription factors are proteins that bind to a specific DNA sequence to control the transcription of specific genes. NAC, WRKY, and MYB transcription factors are central players in modulating transcriptional changes under hormonal control [29, 30]. In our RNA-seq dataset, we identified many NACs, MYBs, and WRKYs from the enriched hormone related GO-terms. As the NACs, MYBs, and WRKYs are the central players involved in the hormonal regulation, we hypothesized that they might play an important role in the HT-induced floral abortion. To quantify the correlation of these central players and the hormone related DEGs underlying floral abortion, we performed PLS-SEM analysis using the RNA-seq data and constructed a reliable model quantifying the relationship among hormone related DEGs, the NACs, MYBs, WRKYs, and the flowering-related DEGs. Interestingly, the model showed that IAA related DEGs were connected with flowering-related DEGs indirectly while JA related DEGs made the contribution directly. Other hormone-related genes might affect it in both of the two pathways. On the whole, the effects of hormone-related DEGs on NACs, MYBs, WRKYs were stronger than those on flowering-related DEGs, suggesting that the hormone related DEGs were more likely to affect flowering-related DEGs indirectly through NACs, MYBs, and WRKYs. Hence NACs, MYBs, WRKYs may act as central players involved in the HT-induced floral abortion underlying hormonal control. Using PLS-SEM analysis, we can figure out which one is more influential between the direct pathway and indirect pathways. Our study shows that it is possible to reveal the relationship between biological processes based on RNA-seq data using PLS-SEM.
From the enriched temperature responsive go term, we identify a homology DREB1A, encoding Dehydration-responsive Element-binding Protein. DREB1A is an APETALA2/ethylene-responsive element-binding factor (AP2/ERF)-type transcription factor. It specifically binds dehydration-responsive elements (DREs) to enhance late embryogenesis-abundant (LEA) protein levels and may accumulate under stress-induced dehydration in plant tissues [31,32,33]. DREB1A overexpression causes late flowering in plants . In the present study, we found that expression of the LcDREB1A increased in the shrinking SP. Interestingly, a morphology comparison study on DPs and SPs showed that SPs might at last shrink to be stunt buds . It seems that dehydration happens during this process. Hence, it is likely that LcDREB1A might be involved in the abortion of the floral buds in relation to dehydration.
To study detail expression patterns of the candidate genes in DPs and SPs, 12 genes identified from the flowering, hormone, and temperature response related DEGs were selected. The expression pattern of the inflorescence in ‘Nuomici’ and ‘Huaizhi’ at the four panicle developmental stages under LT and HT were studied. The results showed that LcLHY, LcMED16, LcSKIP20, LcPOD20 in the two cultivars showed similar expression trends, negatively correlated with flowering. LcAIL6 also had similar expression trends in the two cultivars, but positively correlated to flowering. Hence, it is suggested that increased transcript in LcLHY, LcMED16, LcSKIP20, LcPOD20 and decreased transcript in LcAIL6 might be related to the inhibition of floral development in ‘Nuomici’ and ‘Huaizhi’ under HT conditions. Other candidate genes showed a relationship with flowering in ‘Nuomici’ or ‘Huaizhi’, suggesting that our dataset is applicable for identifying genes involved in the HT-induced floral abortion in litchi.
Previous studies have proved that appropriate LT is needed for floral development in litchi [5, 6]. In accordance with these findings, our experiment performed under controlled environmental conditions also showed that the percentage of flowering terminal shoots of the LT-treated trees was significantly higher than that of the HT-treated trees, and that HT induced floral abortion and resulted in poor flowering. In open fields, warm winter and hot spring frequently happen due to climate change and global warming. Floral development of litchi may be inhibited by HT. In practice, control panicle leaf growth can reduce the HT-induced floral abortion by spraying with the ethylene producer ethephon , ROS inducer methyl viologen dichloride hydrate, and the NO donor sodium nitroprusside . In the present studies, we have identified 127 plant hormone-related DEGs potentially involved in the HT-induced floral abortion. Interestingly, it has been proved that there is cross-talk between hormones and ROS/NO during plant growth and development . Whether ethylene, ROS, and NO have any role in controlling the expression of these hormone-related DEGs and whether they inhibit floral abortion by the same mechanism as the LT does needs further detail investigation.
Six RNA-Seq libraries of DP and SP were constructed, and 3.07–8.97 × 106 clean reads were generated. Digital expression of the DPs with that of the SPs was compared, 1320 up-regulated and 981 down-regulated DEGs were identified. From the enriched GO-term, we identified 54 temperature responsive, 23 hormone homeostasis or biosynthesis related, 137 hormone signal transduction or responsive, and 18 flowering-related DEGs. It is suggested HT might affect flowering-related genes via transcription factors NACs, MYBs, and WRKYs through hormonal control. Increased transcript in LcLHY, LcMED16, LcSKIP20, LcPOD20 and decreased transcript in LcAIL6 might be related to the inhibition of floral development. Our studies provided potential genes for the future molecular breeding of new cultivars that can reduce floral abortion under warm climates, and a novel clue to reveal the relationship of biological processes based on the RNA-seq data using PLS-SEM.
Plant materials and experimental procedures
All the air-layered litchi trees (Litchi chinenesis) were cultivated in the experimental orchard of South China Agricultural University (lat. 23o9’40“N, long. 113o21’18”E). The air-layered seedlings were purchased from a market in Conghua District, Guangzhou. Litchi trees were planted in 30-L pots containing loam, mushroom cinder and coconut chaff (v: v: v, 3:1:1). The potted trees were subjected to LT for floral induction in open fields under winter condition.
For the identification of genes involved in floral abortion, four-year-old air-layered ‘Nuomici’ trees with similar size and phenological stage were selected. When panicle primordia emerged, eight trees were transferred to a growth chamber at 12-h photoperiod with a temperature of 18 °C (low temperature, LT) for 40 d to promote panicle development as controls. Another eight trees were transferred to a growth chamber at 12-h photoperiod with a temperature of 26 °C (high temperature, HT) for 2 weeks to induced abortion of panicles. Flowering conditions of the trees at HT or LT were calculated from eight replicated trees, respectively. Percentage of flowering terminal shoots was calculated as the percentage of the flowering terminal shoots to the total terminal shoots in one tree. The DPs under LT and the SPs under HT as shown in Fig. 1 c and d were collected for microstructure observation, or frozen in liquid nitrogen and stored at − 80 °C for RNA extraction and library construction. Each sample was collected from 2 to 3 trees. Three samples of the SP or DP were used for RNA extraction and library construction.
For the studies of floral development in different cultivars under low and high-temperature conditions, 6-year-old air-layered ‘Nuomici’ and ‘Huaizhi’ litchi trees were selected. Once panicle primordia emerged, five replicated trees of ‘Nuomici’, and five replicated trees of ‘Huaizhi’ were transferred to an LT growth chamber described above. Another five replicated trees of ‘Nuomici’, and five replicated trees of ‘Huaizhi’ were transferred to an HT growth chamber as controls. Percentage of flowering terminal shoots, percentage of leafy panicles, length of the panicles, length of the longest axillary panicles were determined. The inflorescences were collected at four stages according to Yang et al. . These four stages under HT could be defined as primary stage, enlarging stage, shrinking stage, and shrunk stage. Those under LT could be defined as primary stage, enlarging stage, early elongating stage, and late elongating stage as shown in Fig. 7. At stage 1 under HT conditions (HT-S1), the axillary panicle primordia and the rudimentary leaves were at their primary stage of development. The axillary panicle primordia were just visible and the leaflets were conglutinated. At stage 2 (HT-S2), the axillary panicle primordia enlarged, the petiole of the rudimentary leaves began to elongate, and each individual leaflet could be seen. At stages 3 (HT-S3) and 4 (HT-S4), the petioles continued to elongate, and the leaflets began to expand, panicle primordia stopped developing. At stage 1 under LT condition (LT-S1), the axillary panicle primordia were just visible and the leaflets were conglutinated. At stage 2 (LT-S2), the axillary panicle primordia enlarged, the petiole of the rudimentary leaf elongated and each individual leaflet was possible to identify. At stage 3 (LT-S3), the axillary panicle primordia continue to enlarge while the petiole of the rudimentary leaves began to bend. At stage 4 (LT-S4), the axillary panicle primordia continued to enlarge and elongate, while the petiole of the rudimentary leaf continued to curve and could be abscised with a gentle touch. The samples under HT were collected every 8–10 d, while those at HT were collected every 2–3 d according to the developmental procedure. Samples were frozen in liquid nitrogen, stored at − 80 °C for total RNA extraction and gene expression analysis. Each sample of the SP or DP for gene expression analysis was collected from 1 to 2 trees.
Buds were vacuum penetrated, fixed with 4% polyformaldehyde for 4 h. Then 4% ethylenediamine was added to the buffer for tissue softening according to the method of Yang et al. . The tissues were dehydrated in a graded ethanol series from 30 to 100% (v/v), with a final change to xylene, and then embedded in wax and sectioned in 8 μm using microtome (Leica RM2235, Nussloch, Germany). The tissues were stained with hematoxylin eosin and observed under a light microscope (Leica DMLB, Bensheim, Germany).
RNA extraction, library construction, and RNA-sequencing
Total RNA extraction was carried out using Plant Total RNA Isolation kit (Huayueyang, Beijing, China). Oligo-dT beads (Qiagen, Valencia, USA) was used to enrich mRNA. The mRNA was fragmented into short fragments using fragmentation buffer, reversed transcribed into cDNA by random primers. Second-strand cDNA was synthesized using DNA polymerase I, RNase H, dNTPs. After that, the cDNA fragments were purified by Qiaquick PCR Purification Kits (Qiagen, Valencia, USA). The purified cDNA fragments were end repaired, added with poly (A), and ligated to Illumina sequencing adapters. The size selected fragments were amplified and sequenced by Illumina HiSeq™ 2500 at BGI (Shenzhen, China). A 50 single-end (SE) module was used. Six libraries including three biological replicates of DP and three of those of the SP were constructed (Table 1).
RNA-Seq data analysis
Reads were mapped to the litchi transcriptome (http://litchidb.genomics.cn, unpublished) using Bowtie2 (version 2.1.0), with parameter settings followed the manual of eXpress (v1.5.1, https://pachterlab.github.io/eXpress/index.html). The Bowtie2 alignments were then subjected to eXpress for read counts and FPKM calculation of each transcript. All the raw data is available at the NCBI Short Read Archive (SRA) under the accession number PRJNA430479.
Differentially expressed genes (DEGs) between treatments were detected by DESeq2 , based on the ‘eff_counts’ matrices output from eXpress. Significant DEGs were restricted with FDR ≤ 0.01. GO term annotation of litchi transcriptome was obtained from a blastp search against the UniProtKB Swiss-Prot database (http://www.uniprot.org/uniprot/). The associated GO terms of the best hit target were assigned to each litchi transcript. The GO annotation was then used for GO term enrichment analysis of DEGs by the R Bioconductor package GOstats .
Quantitative RT-PCR analysis
First-strand cDNA was synthetized using Reverse Transcriptase M-MLV (RNase H-) system (Takara, Dalian, China) from 1 μg extracted RNA. Quantitative real-time polymerase chain reaction (qRT-PCR) primers F1/R1 (Additional file 8: Table S4) were designed by Primer 6.0 (Premier Biosoft, Palo Alto, USA) and synthesized by Sangon Co. Ltd. (Shanghai, China). The litchi homolog β-actin was used as a reference gene (Additional file 8: Table S4). qRT-PCR was performed according to Lu et al.  on a LightCycler480 real-time PCR machine (Roche, Basel, Switzerland). The transcript quantification of the genes was performed in relation to the reference gene and they were calculated by 2-△△CT method . The analyses were conducted with three biological replicates and three technical replicates.
Data were analyzed by SPSS (version 19.0; IBM Corp., Armonk, NY, USA). The differences among treatment means were evaluated by Duncan’s multiple range test at a 0.05 probability level.
To explore the relationships between DEGs and floral abortion, a hypothetical model according to Chen et al.  was specified and analyzed with PLS-SEM with the support of the SmartPLS 2.0 M3 software . The standardized path coefficient values were generated with the PLS algorithm by Path Weighting Scheme using a bootstrapping method to obtain the significance of path coefficients; the Sign Changes were Individual Changes. The Samples during the calculation of bootstrapping method was 5000, and the Cases of ABA, BR, IAA, CTK, JA, SA related genes’ model was 88, 41, 38, 52, 56, 59 respectively.
Cap-binding Protein Subunit
Agamous-like MADS-box protein
AP2-like ethylene-responsive transcription factor
Homeobox-leucine zipper protein
Zinc finger protein CONSTANS-LIKE 12
Differentially expressed gene
Dehydration-responsive element-binding protein 1A
EARLY FLOWERING 4
Heat stress transcription factor A-2
Probable glutathione S-transferase
Late Elongated Hypocotyls
Mediator of RNA Polymerase II Transcription Subunit
Partial Least Squares Structural Equation Modeling
Reactive oxygen species
Glycine-rich RNA-binding protein
Delta (8)-fatty-acid desaturase 2
Zinc-finger homeodomain protein
Liu WW, Kim HJ, Chen HB, Lu XY, Zhou BY. Identification of MV-generated ROS responsive EST clones in floral buds of Litchi chinensis Sonn. Plant Cell Rep. 2013;32:1361–72.
Menzel CM, Simpson DX. Effect of temperature on growth and flowering of litchi (Litchi chinensis Sonn.) cultivars. J Hortic Sci. 1988;63:349–60.
Chen HB, Huang HB. Low temperature requirements for floral induction in lychee. Acta Hort. 2005;665:195–202.
Yang HF, Lu XY, Chen HB, Wang CC, Zhou BY. Low temperature-induced leaf senescence and the expression of senescence-related genes in the panicles of Litchi chinensis. Biol Plant. 2017;61:315–22.
Zhou B, Li N, Zhang Z, Huang X, Chen H, Hu Z, et al. Hydrogen peroxide and nitric oxide promote reproductive growth in Litchi chinensis. Biol Plant. 2012;56:321–9.
Zhou B, Chen H, Huang X, Li N, Hu Z, Gao Z, et al. Rudimentary leaf abortion with the development of panicle in litchi: changes in ultrastructure, antioxidant enzymes and phytohormones. Sci Hortic. 2008;117:288–96.
Zhou B, Huang X, Chen H, Shu W, Hu Z, Liu W, et al. Effects of ethylene on rudimentary leaf and panicle primordium in litchi: antioxidant enzymes, hydrogen peroxide and nitric oxide. Acta Hortic. 2013;975:247–54.
Yang H, Kim HJ, Chen H, Lu Y, Lu X, Wang C, et al. Reactive oxygen species and nitric oxide induce senescence of rudimentary leaves in Litchi chinensis and the expression profiles of the related genes. Hortic. Res. 2018;5:23.
Lu X, Kim H, Zhong S, Chen H, Hu Z, Zhou B. De novo transcriptome assembly for rudimentary leaves in Litchi chinesis Sonn. and identification of differentially expressed genes in response to reactive oxygen species. BMC Genomics. 2014;15:805–19.
Wang C, Lü P, Zhong S, Chen H, Zhou B. LcMCII-1 is involved in the ROS-dependent senescence of the rudimentary leaves of Litchi chinensis. Plant Cell Rep. 2017;36:89–102.
Saha G, Park JI, Jung HJ, Ahmed NU, Kayum MA, Chung MY, et al. Genome-wide identification and characterization of MADS-box family genes related to organ development and stress resistance in Brassica rapa. BMC Genomics. 2015;16:178.
Li W, Liu X, Lu Y. Transcriptome comparison reveals key candidate genes in response to vernalization of oriental lily. BMC Genomics. 2016;17:664.
Zhou Q, Guo JJ, He CT, Shen C, Huang YY, Chen JX, et al. Comparative transcriptome analysis between low- and high-cadmium-accumulating genotypes of Pakchoi (Brassica chinensis L.) in response to cadmium stress. Environ Sci Technol. 2016;50:6485–94.
Shi R, Wang JP, Lin YC, Li QZ, Sun YH, Chen H, et al. Tissue and cell-type co-expression networks of transcription factors and wood component genes in Populus trichocarpa. Planta. 2017;245:927–38.
Hair JF, Ringle CM, Sarstedt M. PLS-SEM: indeed a silver bullet. JMTP. 2011;19:139–52.
Lu X, Li J, Chen H, Hu J, Liu P, Zhou B. RNA-seq analysis of apical meristem reveals integrative regulatory network of ROS and chilling potentially related to flowering in Litchi chinensis. Sci Rep. 2017;7:10183.
Torti S, Fornara F, Vincent C, Andres F, Nordstrom K, Gobel U, et al. Analysis of the Arabidopsis shoot meristem transcriptome during floral transition identifies distinct regulatory patterns and a leucine-rich repeat protein that promotes flowering. Plant Cell. 2012;24:444–62.
Huang W, Pitorre D, Poretska O, Marizzi C, Winter N, Poppenberger B, et al. ALTEREDMERISTEM PROGRAM1 suppresses ectopic stem cell niche formation in the shoot apical meristem in a largely cytokinin-independent manner. Plant Physiol. 2015;167:1471–86.
Elliott RC, Betzner AS, Huttner E, Oakes MP, Tucker WQ, Gerentes D, et al. AINTEGUMENTA, an APETALA2-like gene of Arabidopsis with pleiotropic roles in ovule development and floral organ growth. Plant Cell. 1996;8:155–68.
Krizek BA, Prost V. Macias a. AINTEGUMENTA promotes petal identity and acts as a negative regulator of AGAMOUS. Plant Cell. 2000;12:1357–66.
Nole-Wilson S, Tranby TL, Krizek BA. AINTEGUMENTA-like (AIL) genes are expressed in young tissues and may specify meristematic or division-competent states. Plant Mol Biol. 2005;57:613–28.
Yamaguchi N, Jeong CW, Nole-Wilson S, Krizek BA, Wagner D. AINTEGUMENTA and AINTEGUMENTA-LIKE6/PLETHORA3 induce LEAFY expression in response to auxin to promote the onset of flower formation in Arabidopsis. Plant Physiol. 2016;170:283–93.
Mizoguchi T, Wheatley K, Hanzawa Y, Wright L, Mizoguchi M, Song HR, et al. LHY and CCA1 are partially redundant genes required to maintain circadian rhythms in Arabidopsis. Dev Cell. 2002;2:62–641.
Fujiwara S, Oda A, Yoshida R, Niinuma K, Miyata K, Tomozoe Y, et al. Circadian clock proteins LHY and CCA1 regulate SVP protein accumulation to control flowering in Arabidopsis. Plant Cell. 2008;20:2960–71.
Janda M, Planchais S, Djafi N, Martinec J, Burketova L, Valentova O, et al. Phosphoglycerolipids are master players in plant hormone signal transduction. Plant Cell Rep. 2013;32:839–51.
Ding X, Cao Y, Huang L, Zhao J, Xu C, Li X, et al. Activation of the indole-3-acetic acid-amido synthetase GH3-8 suppresses expansin expression and promotes salicylate- and jasmonate-independent basal immunity in rice. Plant Cell. 2008;20:228–40.
Choe S, Dilkes BP, Gregory BD, Ross AS, Yuan H, Noguchi T, et al. The Arabidopsis dwarf1 mutant is defective in the conversion of 24-Methylenecholesterol to campesterol in brassinosteroid biosynthesis. Plant Physiol. 1999;119:897–907.
Jirage D, Tootle TL, Reuber TL, Frost LN, Feys BJ, Parker JE, et al. Arabidopsis thaliana PAD4 encodes a lipase-like gene that is important for salicylic acid signaling. PNAS. 1999;96:13583–8.
Caarls L, Pieterse CMJ, Van Wees SCM. How salicylic acid takes transcriptional control over jasmonic acid signaling. Front Plant Sci. 2015;6:170.
Kim HJ, Nam HG, Lim PO. Regulatory network of NAC transcription factors in leaf senescence. Cur Opin Plant Biol. 2016;33:48–56.
Kasuga M, Liu Q, Miura S, Yamaguchi-Shinozaki K, Shinozaki K. Improving plant drought, salt, and freezing tolerance by gene transfer of a single stress-inducible transcription factor. Nat Biotech. 1999;17:287–91.
Maruyama K, Takeda M, Kidokoro S, Yamada K, Sakuma Y, Urano K, et al. Metabolic pathways involved in cold acclimation identified by integrated analysis of metabolites and transcripts regulated by DREB1A and DREB2A. Plant Physiol. 2009;150:1972–80.
Maruyama K, Todaka D, Mizoi J, Yoshida T, Kidokoro S, Matsukura S, et al. Identification of cis-acting promoter elements in cold- and dehydration-induced transcriptional pathways in Arabidopsis, rice, and soybean. DNA Res. 2012;19:37–49.
Kudo M, Kidokoro S, Yoshida T, Mizoi J, Todaka D, Fernie AR, et al. Double overexpression of DREB and PIF transcription factors improves drought stress tolerance and cell elongation in transgenic plants. Plant Biotech J. 2017;15:458–71.
Kocsya G, Tari I, Vankovád R, Zechmanne B, Gulyás Z, Poórc P, et al. Redox control of plant growth and development. Plant Sci. 2013;211:77–91.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11:R106.
Falcon S, Gentleman R. Using GOstats to test gene lists for GO term association. Bioinformatics. 2006;23:257–8.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2-DELTADELTACT method. Methods. 2001;25:402–8.
Chen WL, Li J, Zhu HH, Xu PY, Chen JZ, Yao Q. Arbuscular mycorrhizal fungus enhances lateral root formation in Poncirus trifoliata (L.) as revealed by RNA-Seq analysis. Front Plant Sci. 2017;8:2039.
Ringle CM, Wende S, Will S. Smart PLS 2.0 (M3) Beta. Hamburg; 2005. www.smartpls.de.
This research was supported by the National Natural Science Foundation of China (grant nos. 31772249 and 31572080), the National Litchi and Longan Research System (project no. CARS-33), and the Innovation Team Project of the Department of Education of Guangdong Province (grant no. 2016KCXTD 011). The funding bodies had no roles in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Availability of data and materials
All sequence data have been deposited in NCBI Short Read Archive (SRA) under the accession number PRJNA430479.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. Correlation between qRT-PCR and RNA-seq. Scatter plots represent the fold-changes in the gene expression levels of SP compared to that of DP. (PDF 58 kb)
Figure S2. Correlation analysis of the DPs and the SPs. (PDF 143 kb)
Table S1. Ten top enriched GO terms of the DEGs in biological process (BP), cellular component (CC), and molecular function (MF) ontologies. (PDF 52 kb)
Figure S3. Heat map diagram showing the up-regulated expression profiles of the hormone signaling related differentially expressed genes (DEGs) or hormone responsive DEGs. Litchi trees when panicle primordia emerged were transferred to a growth chamber at 12-h photoperiod with a temperature of 18 °C (LT) to encourage floral development. The other trees were transferred to a growth chamber at 12-h photoperiod with a temperature of 26 °C (HT) to induced floral abortion. FPKM values of the developing panicle (DP) and the shrinking panicle (SP) were normalized to Z-score. (PDF 952 kb)
Figure S4. Heat map diagram showing the down-regulated expression profiles of the hormone signaling related and hormone responsive differentially expressed genes (DEGs). Litchi trees when panicle primordia emerged were transferred to a growth chamber at 12-h photoperiod with a temperature of 18 °C (LT) to encourage panicle development. The other trees were transferred to a growth chamber at 12-h photoperiod with a temperature of 26 °C (HT) to induced floral abortion. FPKM values of the developing panicle (DP) and the shrinking panicle (SP) were normalized to Z-score. (PDF 1712 kb)
Table S2. Total effects and Bootstrapping analysis in PLS-SEM. T statistics higher than 1.96 were significant at 5% according to Hair et al. . (PDF 51 kb)
Table S3. PLS-SEM results quality criteria. All latent variables are significant and goodness-of-fit measures when Average variance extracted (AVE; Indicator for converge validity), and composite reliability (indicator for internal consistency reliability) are equal or higher than 0.5 and 0.7 according to Hair et al. . (PDF 100 kb)
Table S4. Primer sequences of the reference gene and candidate genes for qRT-PCR. (PDF 78 kb)
About this article
Cite this article
Liu, H., Wang, C., Chen, H. et al. Genome-wide transcriptome analysis reveals the molecular mechanism of high temperature-induced floral abortion in Litchi chinensis. BMC Genomics 20, 127 (2019). https://doi.org/10.1186/s12864-019-5493-8