Differential gene expression between the vigorous and dwarf litchi cultivars based on RNA-Seq transcriptome analysis

Litchi (Litchi chinesis Sonn.) is the most economically significant member of Sapindaceae family, especially in sub-tropical regions. However, its tall tree body often brings many inconveniences to production management. In order to modify the tree size or growth for productivity optimization and simplifying management, it is urgent to reveal the dwarf mechanism of litchi for dwarfing rootstocks or cultivar breeding. However, to date, the mechanisms on litchi dwarfism is still poor known. In the present study, transcriptome profiling were performed on L. chinensis cv. ‘Feizixiao’ (FZX, vigorous cultivar) and ‘Ziniangxi’ (ZNX, dwarf cultivar). A total of 55,810 unigenes were obtained, and 9,190 unigenes were differentially expressed between vigorous and dwarf litchi samples. Gene functional enrichment analysis indicated that the differentially expressed unigenes (DEGs) were related to phytohormone metabolism and signal transduction, and energy metabolism pathways. In particular, GA2ox were only up-regulated in ZNX samples, indicating GA might play an important role in regulating huge difference between vigorous and dwarf litchi cultivars. In addition, the 35S::LcGA2ox transgenic tobacco plants were dwarf and had smaller leaves or branches than wild type plants. Our study provided a series of candidate genes to reveal the mechanism of litchi dwarf.


Introduction
In fruit production, tree architecture requires unique horticultural practices, including grafting, pruning, and training [1]. These practices need to be designed to maximize productivity for a minimum of expense. Due to the high cost of labor, especially in developed countries, the modifying of tree size or growth is critical for productivity optimization and simplifying management [2]. Dwarfing rootstocks and/or interstocks have been available and widely used for fruit and nut trees to create orchards with smaller and easier-to-handle trees [3]. For  the use of dwarfing rootstocks has become very common in important temperate fruit trees like apple, pear, peach, and cherry [4][5][6][7]. Unfortunately, in most tropical and subtropical fruit trees (i.e., litchi, longan), dwarfing rootstocks are not commercially available. Furthermore, the mechanism how rootstocks dwarf fruit trees is not clear [8]. Evergreen subtropical crops such as litchi and mango are often hedged or pruned to control the size of the trees [9]. In order to control tree size, plant growth regulators are also applied [10,11], but it will increase the production cost. Genetic engineering offers a promising approach for developing dwarfing fruit trees to minimize negative efforts. Thus far, the progress on genetic manipulation of tree size has been comparatively limited. The barriers include the large size and the long generation times [12]. Quantitative trait locus (QTL) analysis has been performed to aid the tree size breeding [13,14]. However, QTLs have limitations as molecular markers for the early breeding selection [12]. Therefore, the identification and functional analysis of genes associated with tree size is critical for both conventional breeding and genetic engineering. Fruit size is assumed to be controlled by polygenes and the molecular mechanism is not fully understood [15]. Dwarf mutants are effectively used in identification of dwarfism related genes in many fruit trees. Chen et al. [16] obtained a dwarf mutant of 'Williams' variety of banana and elucidated that GA might play a pivotal role in its dwarfism. Recently, the brachytic dwarfism trait (dw) of peach trees was found due to a nonsense mutation in the gibberellic acid (GA) receptor PpeGID1c [1].
GAs play fundamental functions in plant growth and reducing level of active GAs causes the dwarf phenotype in plants [17]. Therefore, the attempts to alter GAs metabolism and/or signaling have been performed to control plant size [18,19]. GA 20-oxidase (GA20ox), GA 3-oxidase (GA3ox), and GA 2-oxidase (GA2ox) catalyzing later reactions are key enzymes controlling GA biosynthesis [20]. These enzymes are encoded by multigene families having different temporal and spatial expression patterns [21]. Down-regulation of GA20ox and GA3ox resulted in decreasing GA levels and displayed a dwarf phenotype. For instance, the suppressed expression of GA20ox gene in apple caused dwarfism, which was restored by the application of GA 3 [22]. By contrast, overexpression of GA2ox genes, which encode enzymes converting active forms of GAs to inactive forms, also produces dwarf plants [23]. Likewise, overexpression of the DELLA genes, which act to repress GA signaling, leads to dwarfism in apple [24]. In addition to GA, brassinosteroids (BRs) have also been linked to dwarfism [12]. Additionally, crosstalk between GA and other phytohormones (i.e. auxin, BRs, ethylene) play essential roles in plant height control [19].
Litchi (Litchi chinesis Sonn.) is the most economically significant member of Sapindaceae family [25]. It has a very long history in China and is famous for its red skin and juicy sweet aril. The litchi tree is medium to large, which can grow up to 10-12 m or even 20 m [26]. Traditionally, litchi are planted with wide spacing with about 70-80 trees per hectare. Such plantings waste land resource in the early years. What is more, there are problems with harvesting, spraying and protection from birds and bats for these large trees [27]. Another problem is, with V-shaped branches, litchi shoots are easily broken off by strong winds [28]. Therefore, conventional high litchi tree architecture costs vast labor and capital on orchard management, and it is extremely necessary for growers to carry out dense and dwarfing planting to reduce production cost. In our previous study, we found litchi cultivars 'Ziniangxi' and 'Ya1' are two of the rare potential dwarf litchi germplasms [29]. And the biological and anatomic characteristics of stem of 'Ziniangxi' are consistent with what Tombesi et al. [30] and Zorić et al. [31] had reported on peach and cherry, respectively. In this case, 'Ziniangxi' litchi might be an excellent germplasm to develop dwarfing and dense plantation on litchi.
In this study, we chose two litchi cultivars 'Feizixiao' and 'Ziniangxi' with different vigorous levels and sequenced the leaves and apical buds by RNA-seq. By comparing anatomical and differential gene expression in vigorous and dwarf cultivars, we hypothesize that genes related to phytohormones pathways and energy metabolism pathways might play important roles in litchi dwarfism.

Plant materials
Litchi chinensis cv. 'Feizixiao' (FZX) and 'Ziniangxi' (ZNX) were planted in the innovation experimental orchard of Hainan Academy of Agricultural Science, Haikou, China. FZX is a vigorous cultivar, with long, sparse, fragile branches, and large, narrow, deep glossy green leaflets. ZNX is a dwarf cultivar, with thin and open spreading branches. Mature leaves and apical buds were collected from the mature trees of FZX and ZNX, respectively. All these samples with three biological replicates were harvested at the same time to avoid the different transcript due to circadian rhythm factors. Trees used in the experiment were not chemically treated or pruned. All samples were immersed in liquid nitrogen and stored at -80˚C for RNA extraction.

Paraffin section microscopy
Paraffin section microscopy was performed following the protocols described by Chen et al. [32]. The paraffin sections were observed using the photomicroscope.

RNA extraction, cDNA library construction, and sequencing
Total RNA was isolated using the Quick RNA Isolation Kit and treated with DNase I (TaKaRa, Japan) to remove genomic DNA contamination. RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). The mRNA enrichment, mRNA fragmentation, second-strand cDNA synthesis, size selection, PCR amplification and sequencing using an Illumina HiSeq (San Diego, CA, USA) were performed at the Novogene Institute (Novogene, Beijing, China).

Data filtering, de novo assembly and functional annotation
Raw data (raw reads) in fastq format were trimmed and filtered using Trimmomatic v0.33 [33]. High-quality reads were obtained by removing reads containing adapters, reads containing poly-N and low-quality reads. At the same time, the Q20, Q30 and GC content were calculated. High-quality reads were aligned to the SSU and LSU rRNA sequences download form silva database [34] using bwa [35]. Next, rRNA reads were removed by a home-made perl script and clean data were obtained. Reads from all libraries were de novo assembled using Trinity [36] and CD-hit [37] into a gene set that served as the reference for subsequent analysis.

Differential expression analysis and functional enrichment analysis
Gene expression levels were calculated based on the length of the gene and reads count mapped to this gene using the FPKM (fragments per kilobase of transcript sequence per millions base pairs) method [38]. Differential expression analysis for each sequenced library was performed using DESeq [39]. The P values were adjusted using the Benjamini & Hochberg method [40]. The corrected P value of 0.05 and abs |log2(Fold change)| of 1 were set as the threshold for significantly differential expression. GO and KEGG pathway enrichment analysis were performed by TBtools (http://cj-chen.github.io/tbtools). REVIGO was applied to visualize the GO enrichment results [41].

Quantitative real-time PCR analysis
Total RNA was isolated as described above and reverse transcribed with oligo (dT) 18 primers using M-MLV reverse transcriptase (Invitrogen, USA). Transcript levels were analyzed using quantitative RT-PCR with the DyNAmo Flash SYBR Green qPCR kit (Thermo, USA) and the CFX96 qPCR System (Bio-Rad, USA). Gene-specific primers were designed using the Primer 5.0 program (PREMIER Biosoft International, Canada) and listed in S1 Table. All reactions were performed in triplicate with three biological replicates. All reactions were normalized using the Ct values corresponding to Lcactin gene (HQ615689). Unigene expression levels were calculated using the 2 ΔΔCT method [42].

Isolation of LcGA2ox genes and functional analysis in tobacco
The full length of LcGA2ox1-3 genes were amplified by PCR with primers (S1 Table) through high fidelity PCR (Prime STARTM HS DNA polymerase, Takara). The PCR products were digested with BamHI and SacI respectively, and fused into the plant binary vector pBI121 digested by the same enzymes. The generated binary vectors were transferred into Agrobacterium tumefaciens strain LBA4404 using the freeze-thaw method. Transformation of Nicotiana tabacum was performed using the leaf disc method as previously described by Chen et al. [43]. Transgenic tobacco plants were confirmed by PCR using genomic DNA. The transgenic (T1) tobacco plants were selected and used for morphological analyses. For morphological analysis, plant height, stem diameter, leaf size, internode number and length were measured.

Statistical analysis
Statistical analyses were performed with SPSS software (SPSS, Chicago, IL). One-way analysis of variance (ANOVA) was used to evaluate the difference on each sample. Heatmap diagram were analyzed using R software with pheatmap methods. Significant correlations between qRT-PCR and transcriptome data were analyzed with SPSS software using Pearson's correlation as the statistical metric. Significant correlations were considered only when an adjusted P value was lower than 0.05.

Anatomical observation
FZX is a vigorous cultivar, with bigger leaf, compound leaf, and longer terminal bud. In contrast, ZNX is a dwarf cultivar, with small leaf, compound leaf, and shorter terminal bud ( Table 1, Fig 1). According to the compound leaf petiole histological studies, huge differences were observed between FZX and ZNX. Compared with ZNX, the pith part and the ray cell is larger in FZX (Fig 1). Therefore, FZX and ZNX were selected to determine the transcriptional changes between vigorous and dwarf litchi cultivars in the present study.

RNA sequencing and transcript assembly
Twelve libraries were generated using the mRNA from four sample groups: FZX-leaves, FZXapical-buds, ZNX-leaves and ZNX-apical-buds, each with three biological replicates. In addition, equal amounts of mRNA from the twelve samples were pooled, generating the library marked 'pool'. These cDNA libraries were then subjected to Illumina deep sequencing. In total, 196,687,468 paired-end clean reads, each 150 bp in length were obtained from the pool ( Table 2). Each library was represented by at least 13.5 million reads, a tag density sufficient for quantitative analysis of gene expression [44]. The reads were mapped to the reference sequence. There were at least 70% total mapped reads for each sample ( Table 2).
All clean reads were de novo assembled into 150,867 transcripts using Trinity [36]. A total of 55,810 unigenes were obtained after redundancy removal using CD-hit [37]. The length of the unigenes ranged from 300-19,032 bp, with N50 of 2,376 bp. There were also 19

Gene annotation and functional classification
All unigenes were annotated by query against various public databases (NR, COG, and KEGG). As a result, 45,740 (81.96% of 55,810) unigenes were matched to one or more of the databases (Fig 2). Of these unigenes, 18,084 (32.40%) were significantly similar to sequences of Citrus sinensis, and 9,026 (16.17%) and 3,899 (6.99%) unigenes showed high similarity to sequences of Citrus clementina and Theobroma cacao, respectively (Fig 3). GO covers three domains: cellular component, molecular function and biological process. A total of 41,852 unigenes (74.99%) could be assigned to at least one GO term (S2 Fig), and detailed information of enriched GO terms was listed in S2 Table. Among these terms, the most representative terms in the biological process category were metabolic process, cellular process and single-organism process. 'Cell', 'cell part' and 'organelle' were the terms that dominated in the cellular component category. 'Catalytic actibity', 'binding' and 'transporter activity' were the most representative terms in the molecular function category. Only one unigene assigned in 'behavior'. All unigenes were aligned against the COG database for functional prediction and classification. In total, 16,599 (29.74% of 55,810) unigenes were assigned appropriate COG clusters, which could be classified into 25 functional categories (S3 Fig). Among them, the largest category was 'General function prediction only' (12.10% of 55,180), followed by 'signal transduction mechanisms' (9.71% of 55,180), 'posttranslational modification, protein turnover, chaperones' (4.29% of 55,180), 'cell wall/membrane/envelope biogenesis' (3.74% of 55,180), and 'Transcription' (3.61% of 55,180).

Differentially expressed genes (DEGs) between vigorous and dwarf litchi samples
The expression patterns of unigenes between the vigorous and dwarf litchi samples were investigated. Pair wise comparison of the samples revealed many DEGs [|log2Ratio|�1, false discovery rate (FDR)�0.001] at the four libraries (FZX-leaves, FZX-apical-buds, ZNX-leaves, ZNXapical-buds) (Fig 4). A total of 9,190 unigenes were found to be significantly differentially expressed in the pair-wise comparisons between any two samples. There were 2,538 DEGs (1,621 down-regulated and 917 up-regulated) between the FZX-leaves and FZX-apical-buds, and most of them were assigned to 'plant hormone signal transduction' KEGG pathways. There were 3,180 DEGs (1,655 down-regulated and 1,525 up-regulated) between ZNX-leaves and ZNX-apical-buds, and most of them were assigned to 'biosynthesis of other secondary metabolites' KEGG pathways. A total of 3,248 DEGs were detected between FZX-apical-buds and ZNX-apical-buds, with 2199 down-regulated and 1019 up-regulated. Regarding KEGG pathways, most of the DEGs were assigned to 'genetic information process' pathway. There were 1,750 DEGs (878 down-regulated and 872 up-regulated) between FZX-leaves and ZNXleaves. Regarding KEGG pathways, most of the DEGs were assigned to 'oxidative phosphorylation' pathways.

DEGs related to phytohormone metabolism pathways
Phytohormone levels have been reported to be closely correlated with the dwarf and development of plants [20,45]. To investigate the relationship between phytohormones and dwarf in litchi, the unigenes related to phytohormone metabolism were analyzed. In the present study, 65 unigenes in the abscisic acid (ABA) metabolism-related pathway, 115 unigenes in the brassinosteroid metabolism-related pathway, 39 unigenes in the ethylene metabolism-related pathway, 35 unigenes in the GA metabolism-related pathway and 97 unigenes in the cytokinins (CTK) metabolism-related pathway, were identified. In this study, 35 unigenes involved in GA biosynthesis pathways were differentially expressed (Fig 5A). Among them, 18 unigenes were up-regulated in the ZNX-leaves, 20 unigenes were up-regulated in the FZX-leaves, 15 unigenes were up-regulated in the ZNX-apticalbuds, and 17 unigenes were up-regulated in the FZX-aptical-buds. In particular, GA3ox (MSTRG.5806, MSTRG.10448 and MSTRG.51345) were only up-regulated in FZX samples and GA2ox (MSTRG.12960, MSTRG.54748 and MSTRG.14947) were only up-regulated in ZNX samples (Fig 5A), indicating GA might play an important role in the huge difference between vigorous and dwarf litchi cultivars.
It is worth noting that the expression of ethylene related unigenes between FZX and ZNX samples had complete opposite trend, 39 unigenes involved in ethylene biosynthesis pathways were differentially expressed (Fig 5B). Almost all these DEGs were up-regulated in ZNX-leaves (16 SAMs, 12 ACO and 11 ACS) and 26 DEGs were up-regulated in ZNX-aptical-buds (9 SAMs, 8 ACO and 9 ACS). On the contrary, there were only 11 DEGs were up-regulated in FZX-leaves (7 SAMs, 3 ACO and 1ACS) and 19 DEGs were up-regulated FZX-aptical-buds (8 SAMs, 5 ACO and 6 ACS). Moreover, the expression of DEGs related to indole-3-acetic acid (IAA), CTK and BR metabolism had similar trend (data not shown). The auxin receptor transport inhibitor response1 (TIR1) and auxin influx carrier protein gene was up-regulated in ZNX-apticalbuds. In addition, the two families of early auxin responsive genes, two out of three GH3 were up-regulated in aptical-buds and one SAUR was only up-regulated in FZX-aptical-buds. Of the DEGs related to CTK signaling, the levels of one type-A response regulator gene increased in ZNX-aptical-buds.

DEGs related to energy metabolism pathways
In this study, a total of 69 DEGs were related to energy metabolism pathways (Fig 6B). Among them, there were 7 ATP synthase subunit unigenes, 5chlorophyll A/B binding protein unigenes, 2 fructose-1,6-bisphosphatase unigenes, 2 fructose-bisphosphate aldolase unigenes, 6 NADH dehydrogenase unigenes, 2 PSII core complex proteins unigenes, 2 PSI reaction center subunit II unigenes, and 2 plasma membrane ATPase-like unigenes. The great majority unigenes were up-regulated in FZX samples, especially in FZX leaves sample (Fig 6B).  (Table 3). There was a high correlation between the qRT-PCR data and the RNA-seq data by linear regression analysis (S5 Fig).

Functional characterization of LcGA2oxs by stable expression in tobacco
As mentioned above, GA2ox (MSTRG.12960, MSTRG.54748 and MSTRG.14947) were only up-regulated in ZNX samples (Fig 5A). To investigate the role of LcGA2ox in litchi dwarfism, their function was investigated by stable expression in tobacco. Three LcGA2ox genes,  .14947), under the expression of the 35S promoter were transformed into tobacco. The existence of introduced LcGA2ox genes was confirmed by PCR. Thirty-two, eleven, and nineteen independent PCR-positive transgenic T1 plants were obtained for LcGA2ox1, LcGA2ox2, and LcGA2ox3, respectively (Fig 7). All transgenic plants were phenotypically distinguishable from wild type plants (Fig 8). Firstly, these 35S::LcGA2ox transgenic plants were smaller than wild type plants (Fig 8A). Secondly, 35S::LcGA2ox transgenic plants flowered later than the wild type plants. When the wild type plants flowered, they were over 80 cm high, whereas 35S:: LcGA2ox transgenic plants were only 6 cm high ( Fig 8B). As many lines had a typical dwarf phenotype, three independent transgenic lines of each 35S::LcGA2ox transgenic plants with a significant change of plant height were selected for further investigation after 6 months of planting. As shown in Table 4, six phenotypes including days to flowering, plant height, stem diameter, leaf area, internode length, and internode number were measured. In particular, the plant height of transgenic tobacco 35S::LcGA2ox3#3 was only 3.28 cm. Besides plant height, leaf area was smaller in in transgenic plants compared with the wildtype (Table 4). These results showed that ectopic overexpression of LcGA2ox genes in tobacco leads to the dwarf phenotype in transgenic plants.

Discussion
Dwarfism is a desirable characteristic for many agricultural plants. Dwarf can reduce lodging and increase harvest index in grain crops. The major factor of the success of the Green Revolution was the dwarf wheat (Triticum aestivum) and rice (Oryza sativa) cultivars breeding [46]. In the modern fruit culture and production, dwarfing and close planting are very important targets. Farmers often use dwarfing rootstocks or dwarf cultivars to develop suitable germplasms and achieve more profits [47]. Many dwarfing apple rootstocks are widely used in production, but these resources are limited in litchi cultivation [48]. Table 3. Confirmation of RNA-Seq expression profiles with qRT-PCR. Lcactin was used as reference gene to normalize gene expression levels under identical conditions.

Unigenes ID
Annotation RNA-Seq (log2FoldChange) qRCR (log2FoldChange) In our previous study, ZNX is identified as a dwarf litchi germplasm [29]. In the present study, the anatomical observation supported this view (Fig 1). Therefore, ZNX might be an excellent germplasm to study the mechanisms of dwarf in litchi. In this study, transcriptome analysis was adopted to understand the global molecular events of differentially expressed genes between vigorous (FZX) and dwarf (ZNX) litchi cultivars. RNA-Seq technology was used to profile the litchi buds and leaves of vigorous and dwarf litchi cultivars transcriptome on the Illumina HiSeqTM 2500 platform, and approximately 197 million paired-end clean reads were obtained. 45,740 out of 55,810 unigenes were successfully annotated against public databases, suggesting their relatively conserved functions. A total of 9,190 unigenes were found to be significantly differentially expressed in the pair-wise comparisons between any two samples (Fig 4), which would provide useful information on the litchi dwarf mechanism.

Anatomical characteristics
Several indicators, such as palisade/spongy ratio, wood/bark ratio, vessel diameter and vessel density, have been used to distinguish dwarfing rootstocks [31,49,50,59,60]. In this study, we also found there were huge differences between vigorous and dwarf litchi cultivars (Table 1, Fig 1). Compared with FZX, the pith part and the ray cell of ZNX is smaller, indicating that ZNX is a potential dwarf litchi cultivar.

Phytohormone related pathways
Among the different factors responsible for plant dwarfism, phytohormone was the main reason. GA and BR have been extensively studied in determining plant height [20,45]. In plants, the flux of bioactive GAs is regulated by the balance between their biosynthesis rates and deactivation. 2-oxidation pathway has been suggested to be a major mechanism for GA inactivation [61][62][63]. Thus, overexpression of GA2ox results in dwarfism [64][65][66]. On the other hand, the suppression of GA2ox expression leads to tall and slender phenotypes [63,67]. In this study, 35 unigenes involved in GA biosynthesis pathways were differentially expressed ( Fig  5A). In particular, GA3ox (MSTRG.5806, MSTRG.10448 and MSTRG.51345) were only upregulated in FZX samples GA2ox (MSTRG.12960, MSTRG.14947 and MSTRG.10600) were only up-regulated in ZNX samples. Additional, ectopic overexpression of LcGA2ox genes in tobacco resulted in a dominant dwarf and later flower phenotype (Fig 8). Our results indicate that the differentially expressed genes involved in GA biosynthesis pathways may result in the difference in growth vigor between FZX and ZNX.
There are few reports about ethylene involving in dwarfism. However, the expression of ethylene related unigenes between FZX and ZNX samples had complete opposite trend in the present study (Fig 5B). Amost all unigenes were up-regulated in ZNX-leaves, including 16 SAMs, 12 ACO and 11 ACS. The correlation between ethylene and litchi dwarfing is worth further research. CTK and auxins play particularly significant role in regulating plant growth and development [67][68][69][70]. In this study, most CTK and IAA related unigenes were up-regulated in the aptical-buds samples. Besides, the up-regulated unigenes in FZX were more than that in ZNX.

Energy metabolism pathways
During the process of photosynthesis, plants convert light energy from the sun into chemical energy stored in molecules. The vigorous plants need stronger energy metabolism to support plant growth and development. Increasing studies suggested that sugars not just as nutrients, but also as signal molecules sensing nutrient status and coordinating plant growth and development accordingly [71][72][73][74][75]. In this study, 69 DEGs related to energy metabolism pathways were found (Fig 6B). Moreover, the great majority DEGs were up-regulated in FZX (the vigorous cultivar) samples, especially in FZX leaves sample, indicating the vigorous FZX need stronger energy metabolism to support vigorous growth.

Conclusions
Transcriptome analysis of differentially expressed genes between FZX (vigorous cultivar) and ZNX (dwarf cultivar) revealed interesting genes that might provide some clues revealing the mechanisms of litchi dwarfism. The expression of genes involved in phytohormone related pathways and energy metabolism pathways exist huge differences between vigorous and dwarf litchi cultivars, especially the GA and ethylene related genes. Ectopic overexpression of LcGA2ox genes leads to the dwarf phenotype in transgenic tobacco. In addition, GA and ethylene might interact with other hormone signaling pathways, such as auxin, ABA, and CTK, forming a complex network to regulate plant growth and development.
Supporting information S1