Abstract

   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 [[36]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 [[37]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 [[38]3]. For example, the use of dwarfing
   rootstocks has become very common in important temperate fruit trees
   like apple, pear, peach, and cherry [[39]4–[40]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 [[41]8].
   Evergreen subtropical crops such as litchi and mango are often hedged
   or pruned to control the size of the trees [[42]9]. In order to control
   tree size, plant growth regulators are also applied [[43]10, [44]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 [[45]12]. Quantitative trait locus (QTL) analysis has
   been performed to aid the tree size breeding [[46]13, [47]14]. However,
   QTLs have limitations as molecular markers for the early breeding
   selection [[48]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 [[49]15]. Dwarf mutants are effectively used in
   identification of dwarfism related genes in many fruit trees. Chen et
   al. [[50]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 [[51]1].

   GAs play fundamental functions in plant growth and reducing level of
   active GAs causes the dwarf phenotype in plants [[52]17]. Therefore,
   the attempts to alter GAs metabolism and/or signaling have been
   performed to control plant size [[53]18, [54]19]. GA 20-oxidase
   (GA20ox), GA 3-oxidase (GA3ox), and GA 2-oxidase (GA2ox) catalyzing
   later reactions are key enzymes controlling GA biosynthesis [[55]20].
   These enzymes are encoded by multigene families having different
   temporal and spatial expression patterns [[56]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]
   [[57]22]. By contrast, overexpression of GA2ox genes, which encode
   enzymes converting active forms of GAs to inactive forms, also produces
   dwarf plants [[58]23]. Likewise, overexpression of the DELLA genes,
   which act to repress GA signaling, leads to dwarfism in apple [[59]24].
   In addition to GA, brassinosteroids (BRs) have also been linked to
   dwarfism [[60]12]. Additionally, crosstalk between GA and other
   phytohormones (i.e. auxin, BRs, ethylene) play essential roles in plant
   height control [[61]19].

   Litchi (Litchi chinesis Sonn.) is the most economically significant
   member of Sapindaceae family [[62]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
   [[63]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
   [[64]27]. Another problem is, with V-shaped branches, litchi shoots are
   easily broken off by strong winds [[65]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 [[66]29]. And the biological and
   anatomic characteristics of stem of ‘Ziniangxi’ are consistent with
   what Tombesi et al. [[67]30] and Zorić et al. [[68]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.

Materials and methods

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. [[69]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 [[70]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 [[71]34] using bwa [[72]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 [[73]36] and CD-hit [[74]37] into a gene set that served as the
   reference for subsequent analysis.

   Gene function was annotated based on the highest similarity in the
   following databases: Nr ([75]http://www.ncbi.nlm.nih.gov, NCBI
   non-redundant protein sequences, with a cut-off e-value of 1e^-5), COG
   ([76]http://www.ncbi.nlm.nih.gov/COG, Clusters of Orthologous Groups of
   proteins), KO ([77]http://www.genome.jp/kegg, KEGG Orthology database).
   GO (Gene Ontology) functional annotation was performed using Blast2GO
   (v2.5.0) software.

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
   [[78]38]. Differential expression analysis for each sequenced library
   was performed using DESeq [[79]39]. The P values were adjusted using
   the Benjamini & Hochberg method [[80]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
   ([81]http://cj-chen.github.io/tbtools). REVIGO was applied to visualize
   the GO enrichment results [[82]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 [83]S1 Table. All reactions were performed in triplicate with three
   biological replicates. All reactions were normalized using the Ct
   values corresponding to Lcactin gene ([84]HQ615689). Unigene expression
   levels were calculated using the 2^ΔΔCT method [[85]42].

Isolation of LcGA2ox genes and functional analysis in tobacco

   The full length of LcGA2ox1-3 genes were amplified by PCR with primers
   ([86]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. [[87]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.

Results

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 ([88]Table 1, [89]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 ([90]Fig 1). Therefore, FZX
   and ZNX were selected to determine the transcriptional changes between
   vigorous and dwarf litchi cultivars in the present study.

Table 1. Comparasion of shoot and leaf length between ZNX and FZX.

   Cultivar Shoot length (cm) Internode length (mm) Leaf length (mm) Leaf
   width (mm)
   FZX 13.21 ± 2.17b 28.22 ± 4.61b 79.08 ± 7.93b 26.58 ± 3.53b
   ZNX 25.60 ± 5.06a 48.72 ± 12.11a 163.45 ± 23.38a 49.09 ± 7.60a
   [91]Open in a new tab

   Data were mean ± SE of three biological replicates. Different lowercase
   letters indicate significant difference at the 0.05 level based on
   Duncan multiple range test.

Fig 1. The comparison of leaf, compound leaf, terminal bud (upper) and
compound leaf petiole longitudinal sections (lower) between FZX and ZNX
litchi cultivars.

   [92]Fig 1
   [93]Open in a new tab

RNA sequencing and transcript assembly

   Twelve libraries were generated using the mRNA from four sample groups:
   FZX-leaves, FZX-apical-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 ([94]Table 2). Each library was represented
   by at least 13.5 million reads, a tag density sufficient for
   quantitative analysis of gene expression [[95]44]. The reads were
   mapped to the reference sequence. There were at least 70% total mapped
   reads for each sample ([96]Table 2).

Table 2. RNA-seq reads in twelve RNA-seq libraries.

   Sample name Number of input reads Average input read length Uniquely
   mapped reads number Uniquely mapped reads% Percent of reads mapped to
   multiple loci
   FZX-apical-bud1 18,443,351 260 14241724 77.22 19.00
   FZX-apical-bud2 23,002,507 255 17229221 74.90 18.85
   FZX-apical-bud3 14,514,169 270 11070898 76.28 18.77
   FZX-leaves1 20,029,366 273 15309985 76.44 19.65
   FZX-leaves2 15,038,810 271 11355168 75.51 20.00
   FZX-leaves3 16,136,701 271 12208298 75.66 19.62
   ZNX-apical-bud1 12,918,190 270 9272820 71.78 21.60
   ZNX-apical-bud2 16,079,578 269 11411651 70.97 22.15
   ZNX-apical-bud3 15,923,038 269 11364189 71.37 21.95
   ZNX-leaves1 15,635,653 270 11126550 71.16 21.90
   ZNX-leaves2 13,732,058 269 9539948 69.47 18.71
   ZNX-leaves3 15,234,047 271 10904961 71.58 21.87
   [97]Open in a new tab

   All clean reads were de novo assembled into 150,867 transcripts using
   Trinity [[98]36]. A total of 55,810 unigenes were obtained after
   redundancy removal using CD-hit [[99]37]. The length of the unigenes
   ranged from 300–19,032 bp, with N50 of 2,376 bp. There were also 19,874
   unigenes (35.61%), 17,584 unigenes (31.51%) in the length range of
   1,000–2,000 bp and 18,352 unigenes (32.88%) with length >2000 bp
   ([100]S1 Fig). Transdecoder from the Trinity package was used for CDS
   prediction. Approximately 48,660 (87.19%) out of all predicted CDS.

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 ([101]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 ([102]Fig 3).

Fig 2. Distribution of annotation results of Unigenes in Nr, GO, KEGG, COG,
Swissprot database.

   [103]Fig 2
   [104]Open in a new tab

Fig 3. Species distribution of best blastx hits to NR database.

   [105]Fig 3
   [106]Open in a new tab

   Cumulative total numbers of unigenes annotated to NR were shown in the
   outermost circle with the dotted line. The size of the area was in
   proportion to the percentage of best blastx hits to the corresponding
   species.

   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 ([107]S2 Fig), and detailed
   information of enriched GO terms was listed in [108]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 ([109]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).

   To further analyze the litchi transcripts, all the unigenes were
   analyzed with respect to the KEGG pathway database. In total, 25,192
   unigenes were assigned to 139 KEGG pathways ([110]S4 Fig). The pathways
   with the most representation among the unique sequences were the
   plant-pathogen interaction (4,504 unigenes), followed by plant hormone
   signal transduction (2,805 unigenes) and starch and sucrose metabolism
   (2,227 unigenes). Furthermore, some important pathways related to plant
   growth, development, morphogenesis and response to stress stimuli,
   including protein processing in endoplasmic reticulum (1,753 unigenes),
   carbon metabolism (1,677 unigenes), ribosome biogenesis in eukaryotes
   (1,639 unigenes), peroxisome(662 unigenes), ABC transporters (615
   unigenes), zeatin biosynthesis (240 unigenes), citrate cycle (TCA
   cycle) (187 unigenes), photosynthesis (157 unigenes) and others, have
   also been annotated successfully.

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,
   ZNX-apical-buds) ([111]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 ZNX-leaves. Regarding KEGG
   pathways, most of the DEGs were assigned to ‘oxidative phosphorylation’
   pathways.

Fig 4. Differential gene expression profiles based on the library between
vigorous and dwarf litchi samples.

   [112]Fig 4
   [113]Open in a new tab

   (A) The numbers of up- and down-regulated genes in comparisons of the
   FZX-apical-buds vs FZX-leaves, FZX-apical-buds vs ZNX-apical-buds,
   FZX-leaves vs ZNX-leaves, and ZNX-apical-buds vs ZNX-leaves litchi
   samples. (B) Venn diagram showing the comparison of differentially
   expressed genes between any two samples.

   Compared FZX-apical-buds vs ZNX-apical-buds with FZX-leaves vs
   ZNX-leaves, there were 1,478 DEGs only expressed in FZX-apical-buds vs
   ZNX-apical-buds, and 2,976 DEGs only expressed in FZX-leaves and
   ZNX-leaves. Moreover, 271 DEGs were both expressed in these two groups
   ([114]S3 Table).

DEGs related to phytohormone metabolism pathways

   Phytohormone levels have been reported to be closely correlated with
   the dwarf and development of plants [[115]20, [116]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 ([117]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-aptical-buds, 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 ([118]Fig 5A), indicating GA
   might play an important role in the huge difference between vigorous
   and dwarf litchi cultivars.

Fig 5. Heat map diagram of relative gene expression levels of DEGs related to
GA (A) and ethylene (B) biosynthesis.

   Fig 5
   [119]Open in a new tab

   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 ([120]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).

DEGs related to phytohormone signal transduction

   A total of 42 DEGs were identified to involving in plant hormone signal
   transduction pathways ([121]Fig 6A). Seven DEGs (MSTRG.16092,
   MSTRG.33657, MSTRG.37884, MSTRG.37950, MSTRG.44318, MSTRG.50722,
   MSTRG.55744) annotated as PROTEIN PHOSPHATASE 2C (PP2C) involved in ABA
   signal transduction were differentially expressed, and six out of the
   seven were up-regulated in FZX samples, and down-regulated in ZNX
   samples. One DEG (MSTRG.20008) annotated as serine/threonine-protein
   kinase SAPK was down-regulated in ZNX-aptical-buds. In the
   auxin-responsive pathway, six unigenes annotated as auxin response
   factors (ARFs) were differentially expressed. Among them, two unigenes
   (MSTRG.15783, MSTRG.34960) were only up-regulated in aptical-buds, two
   unigenes (MSTRG.46032, MSTRG.46092) were only up-regulated in FZX
   samples. The auxin receptor transport inhibitor response1 (TIR1) and
   auxin influx carrier protein gene was up-regulated in ZNX-aptical-buds.
   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.

Fig 6. Heat map diagram of relative gene expression levels of DEGs related to
plant hormone signal transduction pathways (A) and energy metabolism pathways
(B).

   Fig 6
   [122]Open in a new tab

DEGs related to energy metabolism pathways

   In this study, a total of 69 DEGs were related to energy metabolism
   pathways ([123]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 ([124]Fig
   6B).

qRT-PCR validation of RNA-seq-based gene expression

   qRT-PCR was carried out on 12 genes significantly differentially
   expressed as revealed by RNA-seq. They were MSTRG.15647, MSTRG.33946,
   MSTRG.46080, MSTRG.12329, MSTRG.26332, MSTRG.15171, MSTRG.26097,
   MSTRG.21401 and MSTRG.52145. Overall, the results of qRT-PCR were
   consistent with the RNA-seq data ([125]Table 3). There was a high
   correlation between the qRT-PCR data and the RNA-seq data by linear
   regression analysis ([126]S5 Fig).

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)
   FZX-apical-buds vs ZNX-apical-buds FZX-leaves vs ZNX-leaves
   FZX-apical-buds vs ZNX-apical-buds FZX-leaves vs ZNX-leaves
   MSTRG. 15647 alpha-expansin 3 -1.97 0.04 -1.79 1.99
   MSTRG. 33946 phytochrome B -3.02 -0.25 -0.81 -1.25
   MSTRG. 46080 MYB-related protein 308-like -2.74 0.22 -2.18 -3.06
   MSTRG. 12329 Xylem bark cysteine peptidase -2.46 0.51 -1.29 1.40
   MSTRG. 15171 transcription factor AS1 3.89 0.06 2.67 0.18
   MSTRG. 26332 cysteine proteinase 0.05 0.98 0.16 0.25
   MSTRG. 26097 histone H1-like -0.19 1.33 -0.69 -0.18
   MSTRG. 21401 cysteine protease -7.91 -0.75 -0.67 -1.32
   MSTRG. 52145 DELLA protein GAIP-B-like 17.33 1.11 13.34 0.78
   MSTRG.54748 gibberellin 2-oxidase 3.11 -3.00 2.33 -2.37
   MSTRG.12960 gibberellin 2-oxidase 3.63 3.58 3.27 2.89
   MSTRG.14947 gibberellin 2-oxidase 3.46 4.23 3.24 3.89
   [127]Open in a new tab

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 ([128]Fig 5A). To investigate the
   role of LcGA2ox in litchi dwarfism, their function was investigated by
   stable expression in tobacco. Three LcGA2ox genes, LcGA2ox1
   (MSTRG.54748), LcGA2ox2 (MSTRG.12960), and LcGA2ox3 (MSTRG.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
   ([129]Fig 7). All transgenic plants were phenotypically distinguishable
   from wild type plants ([130]Fig 8). Firstly, these 35S::LcGA2ox
   transgenic plants were smaller than wild type plants ([131]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
   ([132]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 [133]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 ([134]Table 4). These
   results showed that ectopic overexpression of LcGA2ox genes in tobacco
   leads to the dwarf phenotype in transgenic plants.

Fig 7. PCR analysis of 35S::LcGA2ox1 transgenic tobacco plants.

   [135]Fig 7
   [136]Open in a new tab

   M, Marker DL2000. P, positive control. wt, wild type.

Fig 8. The phenotypes of transgenic tobacco lines ectopically expressing
LcGA2ox1, LcGA2ox2 and LcGA2ox3.

   [137]Fig 8
   [138]Open in a new tab

Table 4. Characteristics of transgenic lines ectopically expressing LcGA2ox1,
LcGA2ox2 and LcGA2ox3.

   Line No. of plants Days to flowering Plant height (cm) Stem diameter
   (mm) Leaf area(cm^2) Internode length(cm) Internode number
   W. T. 5 79.72 ± 15.34 b 32.41 ± 3.76 a 9.02 ± 0.27 b 207.90 ± 14.65 a
   1.90 ± 0.50 a 15.60 ± 1.15 bc
   35S::LcGA2ox1
   #9 5 106.72 ± 19.91 ab 20.88 ± 1.43b 7.65 ± 0.76 c 132.47 ± 36.26 bc
   1.30 ± 0.35 bc 17.60 ± 1.53 b
   #10 4 104.36 ± 11.99 ab 21.46 ± 3.01b 6.92 ± 0.80 cd 121.09 ± 0.45 bc
   1.43 ± 0.21 b 15.25 ± 2.52 bc
   #11 5 103.33 ± 15.69 ab 21.72 ± 0.84 b 6.18 ± 0.20 d 103.43 ± 4.64 bc
   1.55 ± 0.13 ab 13.30 ± 1.15 c
   35S::LcGA2ox2
   #7 5 109.27 ± 12.93 a 20.96 ± 2.78 b 7.80 ± 0.57 c 141.30 ± 45.89 bc
   1.30 ± 0.21 bc 17.60 ± 1.53 b
   #10 5 112.38 ± 12.97 a 18.75 ± 1.96 bc 8.78 ± 0.72 b 136.07 ± 25.03 bc
   0.85 ± 0.13 cd 21.30 ± 2.08 a
   #11 5 112.17 ± 18.54 a 16.36 ± 2.60 c 9.92 ± 0.31 a 152.42 ± 48.84 bc
   1.27 ± 0.39 bc 13.30 ± 2.52 c
   35S::LcGA2ox3
   #3 5 120.76 ± 14.11 a 3.28 ± 1.30 d 3.93 ± 0.29 e 94.63 ± 38.98 c 0.25
   ± 0.04 e 15.30 ± 2.52 bc
   #14 4 114.59 ± 19.23 a 6.22 ± 0.85 d 7.30 ± 0.43 c 160.99 ± 1.51 ab
   0.34 ± 0.06 e 17.50 ± 1.53 b
   #18 4 117.69 ± 12.42 a 6.54 ± 1.65 d 7.18 ± 0.21 c 119.04 ± 31.01 bc
   0.43 ± 0.04 de 15.25 ± 0.58 bc
   [139]Open in a new tab

   Data were mean ± SE of three biological replicates. Different lowercase
   letters indicate significant difference at the 0.05 level among 8
   sampling sites based on Duncan multiple range test.

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
   [[140]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 [[141]47]. Many dwarfing apple rootstocks are
   widely used in production, but these resources are limited in litchi
   cultivation [[142]48].

   In our previous study, ZNX is identified as a dwarf litchi germplasm
   [[143]29]. In the present study, the anatomical observation supported
   this view ([144]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 ([145]Fig 4), which would
   provide useful information on the litchi dwarf mechanism.

   Several hypotheses have been suggested to explain dwarf mechanisms in
   plant, such as anatomical [[146]49–[147]51], nutritional
   [[148]52–[149]54], and hormonal [[150]55–[151]58].

Anatomical characteristics

   Several indicators, such as palisade/spongy ratio, wood/bark ratio,
   vessel diameter and vessel density, have been used to distinguish
   dwarfing rootstocks [[152]31, [153]49, [154]50, [155]59, [156]60]. In
   this study, we also found there were huge differences between vigorous
   and dwarf litchi cultivars ([157]Table 1, [158]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 [[159]20, [160]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
   [[161]61–[162]63]. Thus, overexpression of GA2ox results in dwarfism
   [[163]64–[164]66]. On the other hand, the suppression of GA2ox
   expression leads to tall and slender phenotypes [[165]63, [166]67]. In
   this study, 35 unigenes involved in GA biosynthesis pathways were
   differentially expressed ([167]Fig 5A). In particular, GA3ox
   (MSTRG.5806, MSTRG.10448 and MSTRG.51345) were only up-regulated 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 ([168]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 ([169]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 [[170]67–[171]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
   [[172]71–[173]75]. In this study, 69 DEGs related to energy metabolism
   pathways were found ([174]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 Table. Primer sequence of genes.

   (DOC)
   [175]Click here for additional data file.^ (37.5KB, doc)
   S2 Table. The list of enriched GO terms.

   (XLS)
   [176]Click here for additional data file.^ (4.1MB, xls)
   S3 Table. The list of DEGs both expressed in FZX-apical-buds vs
   ZNX-apical-buds and FZX-leaves vs ZNX-leaves.

   (XLS)
   [177]Click here for additional data file.^ (139KB, xls)
   S1 Fig. Length distribution of unigenes.

   The x-axis denoted the length range of all groups. The y-axis denoted
   the number of unigenes in each group.

   (TIF)
   [178]Click here for additional data file.^ (15.3KB, tif)
   S2 Fig. Gene ontology classification of unigenes.

   The number of gene GO terms in each functional subcategory was
   presented as the percentage of GO terms for that subcategory out of the
   total GO terms.

   (TIF)
   [179]Click here for additional data file.^ (13.8MB, tif)
   S3 Fig. COG classification of unigenes.

   The y-axis denoted the number of unigenes in each group. The x-axis
   denoted the functional description of each group. Details were shown in
   the right part of the graph.

   (TIF)
   [180]Click here for additional data file.^ (13.9MB, tif)
   S4 Fig. KEGG classification of unigenes.

   The y-axis denoted the number of unigenes in each group. The x-axis
   denoted subclass of KEGG.

   (TIF)
   [181]Click here for additional data file.^ (46.4KB, tif)
   S5 Fig. Correlation between qRT-PCR and data obtained from
   transcriptome analysis.

   The real-time PCR log2 values (x-axis) were plotted against
   colorationstages (y-axis). **indicates a significant difference at p ≤
   0.01.

   (TIF)
   [182]Click here for additional data file.^ (7.6MB, TIF)

Data Availability

   All of the raw reads are available in the NCBI Sequence Read Archive
   database (Accession Number PRJNA503530).

Funding Statement

   China Litchi and Longan Industry Technology Research System (Project
   no. CARS-32-05), National Natural Science Foundation of China
   (31701885), State Key Laboratory for Conservation and Utilization of
   Subtropical Agro-bioresources (SKLCUSA-b201716), YangFan Innovative &
   Entrepreneurial Research Team Project (No. 2014YT02H013), and the
   Guangdong Litchi Breeding, Propagation and Extension Union Construction
   Project. The funders had no role in study design, data collection and
   analysis, decision to publish, or preparation of the manuscript.

References