Comparative analysis of buds transcriptome and identification of two florigen gene AkFTs in Amorphophallus konjac

Leaves and flowers of Amorphophallus konjac do not develop simultaneously thus unique features can be elucidated through study of flowering transformation in A. konjac. In this study, transcriptome libraries of A. konjac leaf buds (LB) and flower buds (FB) were constructed followed by high-throughput sequencing. A total of 68,906 unigenes with an average length of 920 bp were obtained after library assembly. Out of these genes, 24,622 unigenes had annotation information. A total of 6859 differentially expressed genes (DEGs) were identified through differential expression analysis using LB as control. Notably, 2415 DEGs were upregulated whereas 4444 DEGs were downregulated in the two transcriptomes. Go and KEGG analysis showed that the DEGs belonged to 44 functional categories and were implicated in 98 metabolic pathways and 38 DEGs involved in plant hormone signal transduction. Several genes were mined that may be involved in A. konjac flower bud differentiation and flower organ development. Eight DEGs were selected for verification of RNA-seq results using qRT-PCR analysis. Two FLOWERING LOCUS T (FT) genes named AkFT1 and AkFT2 were identified though homologous analysis may be the florigen gene implicated in modulation of A. konjac flowering. These genes were significantly upregulated in flower buds compared with the expression levels on leaf buds. Overexpression of AkFT genes though heterologous expression in Arabidopsis showed that the transgenics flowered at a very early stage relative to wild type plants. These findings indicate that AkFT1 and AkFT2 function as regulation genes in A. konjac flowering development and the two genes may present similar functions during flowering transition.


Materials and methods
Materials. Amorphophallus konjac K.Koch plants were grown in Sangzhi (Zhangjiajie, China) under natural conditions. RNA was extracted from leaf primordium and flower primordium for transcriptome analysis.
Arabidopsis thaliana ecotype Columbia (sustained in our laboratory) was used as the model plant for gene function analysis. Candidate genes were transformed into Arabidopsis through inflorescence infiltration method. Arabidopsis seeds were surface-sterilized, germinated on a plate containing Murashige and Skoog (MS) medium with 2% (w/v) sucrose and 0.75% (w/v) agar supplemented with 30 mg/L basta to select transgenic plants. Plants were transferred into pots containing a mixture of topsoil and vermiculite (3:1). The plants were then grown in a growth chamber at 25 °C under a 16-h light/8-h dark photoperiod.
Construction, sequencing and analysis of cDNA library. Total RNA was extracted from buds using HiPure HP Plant RNA Kit (Magen, China). cDNA libraries were prepared using NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, USA). 1 μg RNA per sample was used as input material for preparation of libraries. The libraries were sequenced using an Illumina Hiseq 2000 platform at Biomarker Technologies Co, LTD (Beijing, China) and paired-end reads of 2 × 100 bp were generated.
Raw data in fastq format were firstly processed using in-house Perl scripts. Clean data were obtained by removing reads containing adapters and low-quality reads from the raw data. Transcriptome assembly was performed using Trinity software with min_kmer_cov set to 2 by default and all other parameters set default.
Differentially expressed genes in the two libraries were identified using DESeq2 R package (https:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ DESeq2. html). p-value corrected by Benjamini-Hochberg method 27 , "p-value < 0.01 and fold change ≥ 2" were used as the criteria for screening differentially expressed genes (DEGs). GO enrichment analysis was achieved by the topGO R package based on the Kolmogorov-Smirnov test. KEGG pathway enrichment was performed by KOBAS 2.0 software with FDR ≤ 0.05. Identification and cloning of FT genes from A. konjac. FT genes were retrieved from unigene annotation of the two cDNA libraries derived from A. konjac buds. Total RNA was extracted from A. konjac using HiPure HP Plant RNA Kit (Magen, China). RNA samples were reverse transcribed to cDNA using ReverTra Ace ® qPCR RT Master Mix (Toyobo, Japan). Full-length coding sequences of A. konjac FT genes were amplified with Golden Star T6 Super PCR Mix (TsingKe, China). The thermocycling conditions were as follows: denaturation of cDNA at 98 °C for 2 min, 32 cycles of 98 °C for 10 s, 60 °C for 10 s, and 72 °C for 5 s. PCR products were cloned into pClone007 Versatile Simple Vector (TsingKe, China) and sequenced at TsingKe Technologies Co, LTD. Phylogenetic analysis was performed using MEGA 7 software 28 . Conserved motifs in protein sequence were identified using MEME tool (http:// meme-suite. org/) with default parameters 29 . Preparation of AkFT1 and AkFT2 overexpression lines. AkFT1 and AkFT2 gene were cloned into pEGAD vector (sustained in our laboratory) at the AgeI and SmaI site using ClonExpress II One Step Cloning Kit (Vazyme, China). Sequencing was performed at TsingKe Technologies Co, LTD to confirm that the genes were cloned successfully. Expression of AkFT1 and AkFT2 gene was modulated by CaMV35S promoter. Arabidopsis thaliana ecotype Columbia plants were transformed using Agrobacterium tumefaciens strain, GV3101 harboring AkFT1 and AkFT2 overexpression vector using the floral dip method 30 .
A total of 8 individual transgenic lines carrying 35S::AkFT1 and 35S::AkFT2 were established, respectively. Among the two transgenic materials, we selected three independent transgenic lines for analysis, and at least 40 plants were selected for observation and statistical analysis.
Real-time quantitative and semi-quantitative RT-PCR analysis. cDNA sequences of A. konjac were obtained as described above. AkEF1-α was used as the internal reference gene 12  Ethics approval and consent to participate. We confirm that the collection of plant material did not involve any endangered or protected plant species and declare that the work reported here is consistent with the IUCN Policy Statement on Research Involving Species at Risk of Extinction and the current laws of China.

Results
Transcriptome analysis. A. konjac is a unique plant with a single petiole and a compound leaf at the top whereby the leaf and flower do not appear concurrently. Therefore, vegetative or reproductive growth of A. konjac is initiated through leaf buds or flower buds (Fig. 1). Two cDNA libraries were constructed from leaf primordium (LB) and flower primordium (FB) of A. konjac to explore gene expression during the development of leaf buds and flower buds, resulting in 48.52 and 47.95 million raw reads in the two cDNA libraries, respectively. After filtering adapters and low-quality reads from the raw data, approximately 47.79 and 47.23 million clean reads were obtained in the two cDNA libraries. Further, all clean reads were de novo assembled using Trinity software 32 and 68,906 unigenes were obtained from the two cDNA libraries. The average size of unigenes was 920 bp, and the N50 length of unigenes was 1403 bp ( Table 2). About 89.11% of unigenes had a length ranging from 300 to 2000 bp.
Unigenes were analyzed using GO, KEGG, COG, KOG, eggNOG, Swissprot, Pfam and nr databases to identify key functions of the genes. A total of 24,622 unigenes were annotated, representing only one third of the total number of unigenes. The number of unigenes annotated in each database was then determined (Table 3). Notably, 24,246 unigenes were annotated in nr database, accounting for 35.19% of the total number of unigenes, whereas only 6739 unigenes had annotation information in COG database. Functional annotation of genes expressed in leaf buds and flower buds further enriched the gene pool of A. konjac.

Differential expression analysis of genes in leaf buds and flower buds.
With LB as the control, 6859 significant DEGs were selected with 2415 DEGs were upregulated whereas 4444 DEGs were downregulated (Fig. 2a). Functional annotation was performed on the identified DEGs according to the expression levels of the genes in two libraries. A total of 2908 DEGs had annotation information (Table 4). Analysis showed that the nr database had the highest number of DEGs with annotated information with 2842 DEGs.
GO functional analysis showed that 1381 DEGs were enriched in 44 classes of three major categories (biological process, cell composition and molecular function) (Fig. 2b). Some DEGs (672) were annotated as metabolic processes which was the most representative class under the "biological process" category. Some DEGs (517) were annotated as "cell part", which was the most significantly enriched term in "cellular component" category. Under "molecular function", DEGs were mainly involved in "binding" (611) and "catalytic activity" (639) processes. TopGO software was used to explore enrichment of DEGs 33 . The top 10 GO terms with the significant enrichment Table 1. Sequences of specific primers.

Name
Forward Reverse EF1-α  AAG TTC CTG AAG AAT GGC GAT  GTC CCT CAC GGC AAA CCT ACC   GA2ox  GTC AAC CTC CAC GCC AAG CAT  CCA GCC TCG ACA TTA CGT CCAG   AkFT1  CGT CAC CAA TGG CTC CGA GTTC  TCC ACC AGC ACG AGT GTG TAGG   AkFT2  TGG ACC CCT TCA CAA GGA CT  GGT CCT GAG ATC GTT GCC TC   WUS  GGT AAT GGC TGT GGT GGC TCTG  CGG CAT TGC TGT TGG CTC CA   EJ2  CAA TCG CCA ACC TGC TCA CTCA  TGG GCA AGT GTT TCT GGG TTCA   GAI1  ATC GGC TCA GCA GCA GCA GTA  TGA TGA GGC GGA GGC AAT GGT   AGL30  TCT TGG TCT      www.nature.com/scientificreports/ of DEGs are presented in Table 5. The top 3 GO terms included DNA integration, structural constituent of ribosome and RNA-dependent DNA biosynthetic process. Moreover, 862 DEGs were annotated into 98 metabolic pathways in KEGG pathway analysis (Fig. 2c). The most significantly enriched pathways were ribosome (80 DEGs), plant hormone signal transduction (38 DEGs) and phenylpropanoid biosynthesis (38 DEGs). Some DEGs (38) were implicated in plant hormone signal transduction and may be involved in flower bud differentiation and floral organ development of A. konjac, and studies should further explore their functions. These metabolic pathways provide a molecular foundation for studying the specific processes involved in leaf bud and flower bud development of A. konjac. KEGG pathway enrichment analysis showed that the top 20 pathways associated with high number of DEGs included ribosome, phenylpropanoid biosynthesis and the pentose and glucuronate interconversions (Fig. 2d). These metabolic pathways are mainly involved in synthesis of organic matter and energy transfer, indicating that there are some differences in material and energy requirements between leaf bud and flower bud development.
DEGs related to flowering. Genes implicated in gibberellin synthesis or flowering signaling pathway were significantly differentially expressed in the two transcriptomes. c75063.graph_c0 and c82483.graph_c0 are homologous to GAI gene and their expression was downregulated in flower buds. The log 2 (fold change) of the two GAI genes was − 6.47 and − 6.70, respectively. c76528.graph_c0 which was homologous to GA20ox gene was highly expressed in flower buds. The log 2 (fold change) of GA20ox gene was 8.26. Expression level of c63309. graph_c0 and c74067.graph_c0 which were homologous to GA2ox was lower in flower buds compared with the expression level in leaf bud. The log2 (fold change) of the two GA2ox genes was − 3.77 and − 6.21, respectively. 4 SPL homologous genes were identified from DEGs and their expression level was significantly high in flower buds relative to the expression level in leaf buds. c73015.graph_c0 and c100034.graph_c0 were homologous FT genes and had high expression level in flower buds. The log 2 (fold change) of the two FT genes was 4.73 and 7.53, respectively. Nine MADS-box transcription factor genes were identified from DEGs, with 7 upregulated genes and 2 downregulated genes in flower buds. Notably, several DEGs which showed specific expression profile in flower buds were implicated in floral development.
Other phytohormonal-related gene expression differences during flowering of A. konjac. Phytohormones participate in a variety of physiological and biochemical processes, and are involved in regulation of growth and development of plants. Expression patterns of phytohormone biosynthesis and signal transduction related genes in leaf buds and flower buds were analyzed to determine the regulatory effect of other phytohormones except GA on the flowering of A. konjac. Expression of auxin biosynthesis gene, YUCCA4; auxin transporter protein gene, AUX1; auxin responsive genes IAA9, IAA10 and IAA13 was upregulated in flower buds. In addition, expression of cytokinin dehydrogenase gene, CKX5 was upregulated in flower buds compared with the expression level in leaf buds. Cytokinin response factor gene, AHP was downregulated in flower buds relative to leaf buds. ABA receptor gene, PYL4 and ABA biosynthesis related genes such as CCD8B and NCED1 were downregulated in flower buds compared with leaf buds. Ethylene (ETH) biosynthesis related genes including ACS1, ACS3 and ACS9 were upregulated. The gene that encodes the rate-limiting enzyme in Brassinosteroid (BR) biosynthesis, DET2 was upregulated. Jasmonates acid (JA) biosynthesis related gene, 4CLL6, and repressor of JA responses including TIFY9 and TIFY10A were downregulated in flower buds relative to the expression level in leaf buds. Differential expression of these genes implied that these phytohormones play complex and different roles in development of leaf and flower buds. Further studies should explore the potential regulatory mechanisms.

Verification of expression of DEGs related to flowering of A. konjac. Eight DEGs implicated in
flowering, of which six genes were upregulated and two genes were downregulated, were used for verification of expression levels obtained from transcriptome analysis. The results of quantitative real time-PCR were in agreement with transcriptome analysis results (Fig. 3), indicating that the transcriptome results were reliable.  (Fig. 4). The results showed that c73015.graph_c0 clustered with ZCN8, and c100034. graph_c0 clustered with Hd3a and VRN3. This indicates that PEBP protein of A. konjac was evolutionary related to PEBP proteins from monocots such as rice, wheat and maize. Conserved motifs in these proteins were analyzed using MEME tool 29 to further explore the function of the two FT homologous proteins (Fig. 4). The findings showed that the proteins shared motifs, indicating that FT protein structure is highly conserved in plants, and they have similar functions. These two FT homologous proteins are implicated in promoting flowering of A. konjac. Notably, c73015.graph_c0 was named AkFT1 and c100034.graph_c0 was named AkFT2 according to the above results.   www.nature.com/scientificreports/ are shown in Supplementary Fig. S1 online. Semi quantitative RT-PCR results showed that AkFT1 and AkFT2 high expression levels in transgenic plants, respectively, but not in WT plants (Fig. 5a,b). Rosette leaves of the two types of transgenic plants was smaller and the vegetative growth period was shorter compared with those of WT plants (Fig. 5c,d). Flowering time of WT plants was 25-28 days, and 10-12 rosette leaves were observed during bolting. The number of rosette leaves ranged from 5 to 6 during the bloom period and the flowering time was about 16 days in the 35S::AkFT1 transgenic plants. Moreover, the 35S::AkFT2 transgenic plants took 18-20 days to bloom, and showed 6-7 rosette leaves at this period (Fig. 5e,f). These results indicated that AkFT1 and AkFT2 may play a role in promoting flowering of A. thaliana, respectively.  www.nature.com/scientificreports/ differentiation and flower organ development 38,39 . FT in Arabidopsis is synthesized in the leaf during the day then it is transported into meristem and functions as a co-transcription factor 19 . Therefore, FT is referred as a florigen. Regulation of the duration of light when corm is in the dormancy period can affect flowering in A. bulbifer 8 . However, corms mainly grow under the soil, and photoperiod may not be the main pathway that affects flowering of A. konjac. Exogenous gibberellin application during corm dormancy accelerates A. muelleri blooming 3,40 , indicating that gibberellin pathway is implicated in regulating Amorphophallus flowering. In addition, the corm age affects its flowering pattern. In the present study, gene expression profile of A. konjac leaf buds and flower buds was explored by transcriptome analysis and several DEGs were implicated in gibberellin synthesis and signaling. We obtained 68,906 unigenes from the leaf bud and flower bud transcriptomes, of which only 24,622 unigenes (35.73%) had functional annotation information. The absence of annotation information for some genes may be due to insufficient genome-wide information, limitations of transcriptome analysis, and incomplete information in functional databases. These genes may have unique functions in konjac and deserve further research. A total of 6859 DEGs were identified by comparing the transcriptome of A. konjac leaf buds and flower buds and 2908 DEGs had functional annotation information. Lack of annotation information of several DEGs can be attributed to lack of reference genome and incomplete database information of A. konjac. This indicates that leaf buds and flower buds in A. konjac exhibit different and unique processes and pathways. Several key genes involved in flowering and the possible regulatory pathways in A. konjac were analyzed. Gibberellin homeostasis in plants is achieved by strict regulation of the activities of "activating enzymes" (GA20ox and GA3ox) and "inactivating enzymes" (GA2ox) 41,42 . GA20ox showed low expression levels in leaf buds, however, its showed high expression levels in flower buds. Two GA2ox genes were significantly downregulated in flower buds compared with the expression level in leaf buds. These results indicated that synthesis and degradation of gibberellin significantly affects flowering of A. konjac. DELLA protein is a negative regulator of gibberellin signaling pathway, and a member of GAI-RGA-and-SCR (GRAS) family 43,44 . A. thaliana expresses five DELLA proteins, including GAI, RGA (REPRESSOR OF ga1-3), RGL1 (RGA-LIKE 1), RGL2 and RGL3 45,46 . Two GAI homologous genes were downregulated or even not expressed in flower buds, which may slow flowering of A. konjac. SPL genes regulate the flowering time of A. thaliana through DELLA-dependent and DELLA-independent pathways. Notably, interaction between SPL genes and miR156 can affect flowering through the age pathway 47,48 . Four SPL homologous genes showed high expression level in flower buds, implying that these genes are implicated in promoting the flowering of A. konjac. The expression product of DEHYDRATION-RESPONSIVE ELEMENT-BINDING PRO-TEIN3 (DREB3) is an AP2/EREBP-type transcription factor and overexpression of DREB3 delays flowering in tobacco 49 . DREB3 was highly expressed in leaf buds, but showed low expression levels in flower buds, indicating that DREB3 is a negative regulator of flowering. Most floral meristem identity genes and floral organ identity genes are members of MADS-box transcription factor family 50,51 . Nine MADS-box genes were identified from DEGs with seven upregulated genes and two downregulated genes in flower buds, implying that different MADSbox genes have different expression patterns and may play different roles in flower development. However, further studies should explore the specific mechanism.

Discussion
In addition to Gibberellin, other phytohormones affect the flowering of A. konjac through complex regulatory mechanisms. Previous studies report that auxin plays a key role in development of inflorescence, flower meristem and flower organs [52][53][54] . Cytokinins modulate initiation and development of reproductive organs 55,56 . ETH is involved in regulation of the flowering time of plants 57 . In the current study, significant differences were observed in the expression profiles of several important genes related to auxin, cytokinins, ABA, ETH, JA and Br biosynthesis and signal transduction in leaf and flower buds. Notably, YUCCA4 (Auxin biosynthesis), CCD8B (ABA biosynthesis), NCED1 (ABA biosynthesis), ACS1 (ETH biosynthesis), ACS3 (ETH biosynthesis), ACS9 (ETH biosynthesis), DET2 (BR biosynthesis), 4CLL6 (JA biosynthesis) were differentially expressed in leaf and flower buds. Differential expression of genes implicated in plant hormone signal transduction affects expression of responsive genes, which may ultimately affect the flowering of A. konjac.
The FT gene encodes Phosphatidyl Ethanolamine-binding Protein (PEBP) protein which is a plant flowering integration factor 58,59 . Currently, FT genes from different plants have been reported 21,22,25,[60][61][62] . Some plants have several types of FT genes, and their functions are different. Previous studies report that FT can integrate signals from various flowering pathways to promote plant flowering [63][64][65][66] . FT-like and TERMINAL FLOWER 1(TFL1)-like genes affect several physiological processes in plants, such as seed development and germination in A. thaliana 23,67 , corm formation of potato 68 and bulb development of onion 24 . In this study, two FT homologous genes were isolated from flower buds of A. konjac and were named AkFT1 and AkFT2. Sequence alignment showed that the two FT proteins are high homologous to several plant FT proteins and belonged to FT-like protein. Phylogenetic analysis showed that these genes were highly related to FT protein expressed in monocotyledons, such as rice, wheat and maize. The two FT homologous proteins shared conserved motifs with several FT proteins, indicating that they may have similar biological functions. AkFT1 and AkFT2 were significantly upregulated in flower buds relative to the expression level in leaf buds. This finding was further confirmed by qRT-PCR analysis, indicating that they may be positive regulators of flowering in A. konjac. Significantly, FT protein is transported via the phloem for a long distance, and finally lead to the formation of flowers in the shoot apex 17,69 . This unique feature of A. konjac leaf and flower not developing simultaneously, and the detection of AkFT1 and AkFT2 mRNAs in flower buds, suggests that the expression pattern of FT in Konjac differs from that of other plants. Overexpression of AkFT1 and AkFT2 reduced the vegetative growth period of A. thaliana and accelerated flowering compared with the wild type plants. AkFT1 and AkFT2 genes may play a critical role similar to the function of florigen during flowering transition of A. konjac, however, the precise molecular mechanism should be explored further.
In summary, comprehensive gene expression information of leaf buds and flower buds of A. konjac was obtained through transcriptome analysis. These results showed that some genes are differentially expressed during the development of leaf buds and flower buds. Two FT homologous genes (AkFT1 and AkFT2) were www.nature.com/scientificreports/ identified, which exhibited high expression level in flower buds relative to the expression level in leaf buds. Overexpression of AkFT1 and AkFT2 significantly decreased flowering time of transgenic A. thaliana relative to the wild type plants.