Transcriptome Analysis of Storage Roots and Fibrous Roots of the Traditional Medicinal Herb Callerya speciosa (Champ.) ScHot

Callerya speciosa (Champ.) ScHot is a woody perennial plant in Fabaceae, the roots of which are used medicinally. The storage roots of C. speciosa are derived from fibrous roots, but not all fibrous roots can develop into storage roots. To detect key genes involved in storage roots formation, we performed Illumina sequencing of the C. speciosa storage roots and fibrous roots. De novo assembly resulted in 161,926 unigenes, which were subsequently annotated by BLAST, GO and KEGG analyses. After expression profiling, 4538 differentially expressed genes were identified. The KEGG pathway enrichment analysis revealed changes in the biosynthesis of cytokinin, phenylpropanoid, starch, sucrose, flavone and other secondary metabolites. Transcription factor-related differentially expressed genes (DEGs) were also identified, including such gene families as GRAS, COL, MIKC, ERF, LBD, and NAC. The DEGs related to light signaling, starch, sugar, photohormones and cell wall-loosening might be involved in the formation of storage roots. This study provides the first transcriptome profiling of C. speciosa roots, data that will facilitate future research of root development and metabolites with medicinal value as well as the breeding of C. speciosa.


Introduction
Callerya speciosa (Champ.) ScHot is a woody perennial plant in Fabaceae and distributed over South China and Southeast Asia. Its roots have long been used for food and as a traditional medicinal herb, with properties of toxin removal, heat clearance from the lungs to relieve cough, liver purging and kidney invigoration. More recently, C. speciosa has been widely planted because the wild plants are on the brink of extinction, and rapid propagation has been applied for the manufacture-scale production of seedlings. In previous studies of C. speciosa, we have developed SSR markers to analyze population differences and diversity among natural populations, screened superior variety, and investigated the application of cryopreservation [1,2]. The thickness and starch content are the main factors of quality of the storage roots. The storage roots (SRs) are string-or spindle-shaped and may be lignified and stop expanding under uncertain conditions. The FRs of C. speciosa possess the potential to form SRs, but limited FRs can transfer into SRs, the mechanism by which this occurs is poorly understood.
During evolution, some plants have acquired the ability to differentiate leaves, stems, or roots into storage organs in response to drought or freezing conditions [3]. The mechanisms of storage-organ formation have been investigated in root and tuber crops such as potato (Solanum tuberosum), sweet potato (Ipomoea batatas), cassava and yam. In particular, the tuberization mechanism has been most studied in potato. Homologs of the flowering locus T (FT) family are described as key components of signals inducing the formation of potato tuber and onion bulb [4,5]. In potato, CONSTANS (CO) homologs repress the expression of mobile tuberizing signals under non-inductive long-day conditions (LD). StCOL2, which is expressed at higher levels in potato than other CO homologs, may repress tuberization by directly regulating the transcription of StSP5G, a putative tuberization repressor that may inhibits the expression of the mobile tuberizing signal StSP6A [4,6]. In addition, potato homologs of CYCLING DOF FACTOR 1 (StCDF1), GIGANTEA (GI) and Flavin-binding kelch repeat F-box protein 1 (FKF1) form a complex that degrades StCDF1, thereby repressing the expression StCOL1 and StCOL2 and inducing tuber formation [7,8]. Several plant hormones have been proven to be involved in tuberization in potato [9]. For instance, a decrease in the level of gibberellic acid (GA) is required for tuberization onset, indicating that GA plays an inhibitory role in tuberization [10,11], and given the dramatic change in StARF and StPIN family gene expression during tuberization, it has been suggested that auxin functions as a promoter of tuber formation [12,13]. By promoting cell division during tuberization onset and sink formation, cytokinins may also serve as universal regulators of storage-organ formation [14]. Some other mobile signals have been found to regulate tuberization in potato. miR156 and miR172, two major components of the flowering age pathway, were shown to travel to the stolon via the phloem and modulate tuber formation, with miR156 accumulating in stolons under LD and miR172 during tuberization onset [15,16]. Overexpression of StBEL5 and POTATO HOMEOBOX1 (POTH1) in potato resulted in enhanced tuberization, and movement of the transcripts was demonstrated [17][18][19]. Furthermore, many tuberization-related transcription factors (TFs) such as NAC family genes, MADS-box family genes and knotted-like homeobox genes have also been found in sweet potato [20][21][22][23][24][25].
Several studies have applied transcriptome analysis to detect key genes involved in storageorgan formation of root and tuber crops. Firon et al. [26] performed transcriptome profiling of sweet potato roots, reporting up-regulated starch biosynthesis and down-regulated lignin biosynthesis in SRs compared to FRs. Sun et al. [27] studied tuberous root development in Rehmannia glutinosa using dynamic transcriptome profiling, and Shan et al. [28] used digital gene expression (DGE) tag profiling analysis to identify putative genes involved in photoperiodic tuberization in potato. Additionally, Sojikul et al. [29] revealed phytohormone action during cassava storage root initiation through genome-wide microarray analysis, and Yang et al. [30] employed expression profiling with microarray to demonstrate an active process of glycolysis/gluconeogenesis in cassava storage roots. As there is limited genomic and transcriptomic information for C. speciosa to date, in the present study, we applied Illumina RNA-seq to perform the first transcriptome profiling of C. speciosa roots to find candidate genes involved in SR formation and to provide sequence resources for further analysis.

Sequencing and de novo assembly
To characterize transcriptome differences between SRs and FRs, total RNA was extracted in replicates to prepare cDNA libraries and then subjected to sequencing (S1 Fig). A total of 315,665,026 clean reads of 90 nt were generated (Table 1) and subsequently de novo assembled using Trinity, resulting in 161,926 unigenes with an N50 value of 2107 nt and a mean length of 1285 nt (Table 2). Using the CEGMA pipeline, we checked the completeness of the assembly by similarity searches of 248 conserved eukaryotic core genes [31]. The result indicated that 95.97% of the core genes were completely assembled and 3.81% were partially assembled, suggesting the completeness of the assembly.

Functional annotation of the unigenes
For functional annotations, the assembly unigenes were searched against the NCBI non-redundant protein (Nr), Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), Cluster of Orthologous Groups (COG), NCBI Nucleotide (NT) and SwissProt databases using BLAST with an e-value cutoff of 1e-5. Of the 161,926 assembly unigenes, 125,967 matched to sequences in the databases. According to the BLAST results in the order of Nr, SwissProt, KEGG and COG, we predicted the coding sequences (CDSs) of 109,421 of the unigenes matching to sequences in the protein databases. Of those unigenes for which no sequence was matched, we predicted the CDSs of 8,159 using ESTScan. There are still 44,346 (27.39%) unigenes that have no CDS, some of which may function as non-coding RNAs. Based on the results of BLAST searches against the Nr database, 111,706 unigenes match to sequences of Glycine max (61,295), Medicago truncatula (26,330), Lotus corniculatus var. japonicus (4649), Vitis vinifera (2662), Amygdalus persica (1535), Ricinus communis (905) and other species (14,330) (Fig 1B). Of the 111,706 unigenes, 65,639 (58.76%) showed an e-value less than le-45, and 53,610 (47.99%) showed sequence similarity greater than 80% (Fig 1A). To further characterize the function of unigenes, we applied GO functional classification to unigenes based on the results of BLAST against the Nr database. In total, 87,124 unigenes were classified into 56 functional groups according to the categories of biological process, cellular component and molecular function (Fig 2). Among the 22 functional groups in biological process, cellular process (54,563) and metabolic process (53,339) are the most highly represented groups, followed by single-organism process (36,337), response to stimulus (25,626), biological regulation (21,305), regulation of biological process (19,606), multicellular organismal process (14,433), localization (14,424), and developmental process (14,305), among others. Similar to biological process, unigenes were grouped into 16 sub-groups of molecular function, among  which binding (45,752), catalytic activity (45,038), transporter activity (5818), nucleic acid binding transcription factor activity (2299), and structural molecule activity were the top 5 represented groups.
Based on a BLAST search against the KEGG database, 69,375 unigenes were assigned to KEGG pathways. Metabolic pathways and biosynthesis of secondary metabolites were the most represented pathways, followed by plant-pathogen interaction, plant hormone signal transduction, RNA transport, and spliceosome. Among metabolic sub-pathways, global map was the most represented, followed by carbohydrate metabolism, lipid metabolism, amino acid metabolism, nucleotide metabolism, and biosynthesis of other secondary metabolites ( Fig  3C). For carbohydrate metabolism, sub-pathways of starch and sucrose metabolism, pyruvate metabolism and glycolysis/gluconeogenesis were the most highly represented ( Fig 3A); phenylpropanoid biosynthesis, flavonoid biosynthesis, stilbenoid, diarylheptanoid and gingerol biosynthesis, flavone and flavonol biosynthesis and isoflavonoid biosynthesis were the top five represented sub-pathways of metabolism of other secondary metabolites ( Fig 3B). Otherwise, some researches indicated that that many likely medical substances such as isoliquiritigenin, pterocapan, daucosterol, beta-Sitosterol, stigmasterol, chalcone, flavonoid, polysaccharide existed in roots of C. speciose, of which genes and pathways were identified in KEGG annotations [32,33].

Differential expression analysis
We first estimated the expression level of unigenes using fragments per kilobase of transcript per million fragments mapped (FPKM) values. The SR libraries and FR libraries contained 142,510, and 156,700 unigenes, with 19,416 specifically expressed in the SR libraries and 5,266 in the FR libraries. Using the NOIseq method proposed by Tarazona et al. [34] for differential expression analysis with |log 2 (SR/FR)|1 and probability0.8, we identified 4538 differentially expressed genes, among which 491 were up-regulated and 4047 down-regulated (S1 Table).
To explore the function of DEGs, 1579 of 4538 DEGs were assigned KEGG annotations and then subjected to KEGG pathway enrichment analysis. As a result, DEGs were classified to 118 pathways, 46 of which showed significant enrichment.

Identification of transcription factor genes among DEGs
To identify transcription factor (TF)-associated unigenes, BLASTx was used to search PlantTFDB for DEG sequences [35], with 1597 DEGs classified into 54 transcription factor gene families. The ERF (117) gene family was the most represented, followed by WRKY (111)

Identification of DEGs involved in formation of storage roots
We identified DEGs might be involved in SR formation (Table 3, S2 Table). First, we examined unigenes might be involved in flowering timing, revealing homologs of CO, FT, TERMINAL FLOWER 1 (TFL1), and Flowering-promoting factor 1 (FDF1). Homologs of FT and TFL1, CONSTANS-LIKE 2 and CONSTANS-LIKE 5 were up-regulated, whereas homologs of FDF1 and EARLY FLOWERING 4 (ELF3) were down-regulated. Moreover, many TFs might be involved in flower and root development were identified. A large number of plant-specific transcription factors belong to NAC domain-containing protein families, some of which are involved in cell division and expansion, such as NAC21/22, NAC29 and NAC043; these homologs were down-regulated in SRs. LOB-domain (LBD) genes comprise a large family of plantspecific and DNA-binding transcription factors, the functions of which are largely unknown [36], and homologs of LBD4 and LBD5 were up-regulated but homologs of LBD13, LBD38 and LBD41 down-regulated. RAP1 belongs to the BHLH family of transcription factors, which are involved in anthocyanin biosynthesis and phytohormone signaling, and the homolog was found to be up-regulated in SRs. ARFs are transcription factors that response to auxin and function in promoting flowering, stamen development, and floral organ abscission as well as repressing cell division and organ growth [37]; homologs of ARF4, ARF5, ARF6, ARF9 and ARF18 were down-regulated. In the meanwhile, the homologs of auxin influx transporter LAX1, LAX2, LAX3 and LAX5 were down-regulated in SRs.

qRT-PCR validation of DEGs
Fifteen DEGs were selected for qRT-PCR, all of which were found to be up-regulated or down-regulated in agreement with the differential expression analysis by RNA-seq (Fig 6).  CL11331.Contig2_All, CL13177.Contig2_All, CL14063.Contig6_All and Unigene25102_All were up-regulated more than 25-fold in SRs, and the relative was up-regulated expression level less than 7-fold; the qRT-PCR-based results for Unigene41230_All, Unigene27068_All and Unigene18836_All also revealed a lower degree of relative down-regulation. This result

Discussion
This study provides a comparative transcriptome analysis between C. speciosa SRs and FRs. A total of 4538 DEGs were identified. KEGG pathway enrichment analysis revealed that DEGs were enriched in pathways such as starch and sucrose metabolism, plant hormone signal transduction, and phenylpropanoid biosynthesis, indicating changes in carbon flow and phytohormone signaling might happened during the formation of SRs of C. speciosa. DEGs were also enriched into the pathways Flavonoid biosynthesis, Flavone and flavonol biosynthesis, Isoflavonoid biosynthesis, Benzoxazinoid biosynthesis, Isoquinoline alkaloid biosynthesis, Indole alkaloid biosynthesis, and Glucosinolate biosynthesis, demonstrating the different secondary metabolite biosynthesis between SRs and FRs. Furthermore, DEGs related to transcription

Unigenes involved in light signaling and flowering
Appropriate photoperiod is acquired for initiation of storage organ in potato, lotus and onion [38][39][40]. And light is necessary for initiation of tuberous roots of R. glutinosa, with the intonation of tuberous roots were completely blocked when shaded with sunshade net [41]. Light signaling related genes COL2, COL5, Early flowering 4 (ELF4) and serine/threonine protein kinase WNK8 differently expressed between SRs and FRs. COL5 affects the expression FT and induce flowering in short-days [42], and COL2 responds to light stimulus but not influence flowering in A. thaliana [43]. ELF4 and WNK8 control flowering time in A. thaliana, with ELF4 involved in circadian rhythms and WNK8 modulating photoperiod response [44,45]. FT and TFL1 both function in flowering time control and tuberization in potato, were up-regulated in SRs of C. speciose [4,46]. Therefore, light signaling might be involved in the formation of SRs of C. speciose, and some flowering related genes expressed in roots may also have distinct functions in root development. Although the storage organs of C. speciose, R. glutinosa, potato, lotus and onion derive from different organs, similar regulation mechanisms may exist, on which light signaling might play a key role. Large populations of mRNAs were found to be able to move between shoots and roots in Arabidapsis thaliana [47]. The mobility of mRNAs related to transcript abundance and halflife [48]. The mRNA of FT is involved in flowering induction by long distance movement in A. thaliana [49]. Some researches indicate that the expression of FT in underground stolons is activated by mobile FT protein in potato [4], which has not been found in other plants. However, there is no evidence that FT mRNA can transport from leaves into roots, although many mRNAs may able to move between shoots and roots in C. speciose. Therefore, further research is needed to determine whether the mRNA of FT in roots of C. speciose is transported from leaves or expresses locally.

Unigenes related to sucrose and starch
A surplus of sugar in plants may be one of the main causes of storage root formation. Indeed, sugars can function as signaling molecular and gene regulators, regulating processes such as growth, flowering transition and tuber formation [50]. Sucrose has been found to promote tuber formation in potato, yam and R. glutinosa [51][52][53]. Sucrose synthase (SUS) converts sucrose into fructose and UDP-glucose, functioning as a sucrose-clearing enzyme [54]. Conversely, Sucrose-phosphate synthase (SPS), the key enzyme catalyzing the synthesis of sucrose, drives sucrose accumulation. The sucrose synthase alleles up-regulated in overall, for that homolog of SUS4 showed an expression level 10-fold than homolog of SUS2, whereas SUS4 upregulated and SUS2 down-regulated in SRs. The up-regulated sucrose synthase may also play a key role in the formation of SR of C. speciose.
The main form of carbohydrate stored by SRs of C. speciose is starch. The starch content of dried SRs is 49.3%, while dry matter content of the SRs is about 67% [55]. The storage organ of C. speciose and sweet potato both are modified roots that accumulated starch. The homologs of ADP-glucose pyrophosphorylase and granule-bound starch synthase, the key enzyme for starch synthesis, and sucrose synthase functioned in both sucrose synthesis and cleavage, up-regulated during SRs formation both in sweet potato and C. speciose, with similar expression pattern during SRs formation [26,56]. The homolog of beta-fructofuranosidase, which functions as sucrose cleavage enzyme, showed a very low expression level, with a mean FPKM value of 0.16 and 0.07 in SRs and FRs of C. speciose respectively. The expression level of beta-fructofuranosidase of sweet potao were very low and decreased during SRs development [26,56], whereas the expression level was higher in SRs than FRs in C. speciose, indicating that beta-fructofuranosidase may play different roles during SRs development between C. speciose and sweet potato. The transcriptome data also suggested that sucrose synthase was predominant sucrose cleavage enzyme for starch accumulation during SRs development of C. speciose, which is similar with sweet potato. Furtherly, study on proteomic of cassava during formation of tuberous roots indicated that 14-3-3 protein up-regulated and play important roles in starch accumulation [57]. But homologs of 14-3-3 like protein C and 14-3-3 protein 7 down-regulated in SRs, showing a possible and different role of 14-3-3 protein during SRs formation and starch accumulation.
Laccases, which are involved in the lignin catabolic process, are required for secondary cell wall lignification [58,59], and Cinnamyl alcohol dehydrogenase (CAD) also participates in lignin biosynthesis by catalyzing the final step of lignin monomer production. Homologs of laccases and CAD were down-regulated, indicating that lignin biosynthesis and lignification might decreased during SRs formation. The homolog of Anthocyanin regulatory C1, the key trans-acting factor required for anthocyanin biosynthesis, and the homolog of Root-specific chalcone synthase were both up-regulated in SRs, suggesting that changes in sugar metabolism may accompany changes in anthocyanin biosynthesis.

Crosstalk between phytohormones, sugars and calcium signaling
Phytohormones play key roles in tuber or tuberous root formation. In potato, GA inhibits tuber formation and promotes stolon elongation [60]; in sweet potato, high cytokinin and auxin levels promote storage root initiation and growth [61,62], and ABA levels positively correlate with the thickening of sweet potato storage roots [63]. GA-, ABA-, ethylene-, auxinrelated biosynthesis unigenes were down-regulated, whereas cytokinin-related unigenes were up-regulated, indicating that cytokinin may participle in the formation of SRs in C. speciosa.
Plants possess systems that link carbon status and development, in wich the Hexokinase (HXK) play a key role. HXK-based signaling interacts negatively with auxin but positively with cytokinin [64,65]. And HXK1 signaling involves extensive crosstalk with plant hormone signaling via interaction with F-actin [66]. The up-regulation of homologs of HXK1, suggested that sugars may affect the formation of SRs partly via the regulation of complex networks of sugar-sensing systems in C. speciosa.
Sugar-based signaling pathways also interact with the calcium signaling pathway. CIPKs function together with calcium-sensing CBLs, comprising a CPIK-CBL network involved in plant signal transduction in response to abiotic stress [67]. AtCIPK4 was found to be induced by salt stress, whereas OsPK4, the rice CIPK4 homolog, was induced by illumination, cytokinin treatment and nutrient deprivation [68,69]. AtCIPK7 is induced by cold and sugars and may function in sugar metabolism through the phosphorylation of sucrose synthase [67,69]. In our study, homologs of CIPK4 and CIPK7 were up-regulated in SRs, indicating that CIPKs may also be involved in regulating the formation of SRs response environment and endogenous cues.

Unigenes involved in cell-loosening and root development
Growth in cell volume begins with selective cell wall-loosening, in which cell wall stress is relaxed, resulting in consequential water uptake and cell enlargement [70]. The growing plant cell wall consists of one or more layers of cellulose microfibrils, with each microfibril directly contacting another or interacting with matrix polymers embedded in a hydrophilic matrix [71]. Expansins meditate acid-induced growth through non-lytic cell wall loosing. Based on sequence based phylogeny, plant expansins can be classed into two large families, EXPA and EXPB, and two small families, EXPLA and EXPLB [71]. In C. speciose, homologs of EXPA1, EXPA10 and EXPA15 are up-regulated and homologs of EXPA6 and EXPLB1 are down-regulated, of which the function need to character in the further research. Several transcription factors have been found to regulate expansins genes. Arabidopsis thaliana homeobox 12 (ATHB12) increases the expression of AtEXP10, promoting cell expansion and leaf growth [72]. In SRs of C. speciose, ATHB12 is down-regulated and the differential expression of expansins may attributable to the homologs of ATHB12.
It has been proposed that the root architecture influences the formation of tuberous root, with the development of lateral roots preventing lignification of the stele in the main axis of the adventitious root, which described in the swelling of sweet potato [26,73]. LBD family genes are involved in root formation and are up-regulated in sweet potato tuberous roots [26,74]. In addition, LBD3 and LBD4 are induced by cytokinin and required for a flat and asymmetrical leaf lamina in Arabidopsis [36]. LBD4 may also be involved in SRs formation of C. speciose, with up-regulation of homologs of LBD4 in SRs. NAC21/NAC22 may function as a transcription activator mediating auxin signaling to promote lateral root development [75], and NAC29 may function in the transition between active cell division and cell expansion [76]. NAC043, also called NAC SECONDERY CELL WALL THICKENING PROMOTING FACTOR 1 (NST1), functions together with NST2 and NST3 in cell wall thickening and cell wall lignification [77][78][79]. The differences of root architecture between C. speciose SRs and FRs might partly resulted from the down-regulation of homologs of NAC21/ NAC22, NAC29 and NAC043 in SRs, which is need further research.

Conclusion
We performed comparative analysis of C. speciose SRs and FRs, discovering a series of genes associates with light signaling, phytohormones, starch, sugar and cell wall-loosening differentially expressed between SRs and FRs, which may be involved in the formation of SRs. However, further research is required to determine the DEGs responsible for SR formation, which will shed lights on the mechanism of SR formation and help in the breeding of C. speciose.

Plant material
Three in vitro-grown seedlings transplanted in a greenhouse for 5 months were selected for the present study. The experimental site is the Institute of Tropical Crop Genetic Resources, Chinese Academy of Tropical Agricultural Sciences (CATAS). The greenhouse are on the conditions of temperature of 20-35°C, relative humidities of 85%-95%, transmittance of 20% and natural photoperiod. The three plants were planted in the nursery bed and watered and fertilized weekly. The potting soil was a mixture of equal parts coconut fibre, sand and turfy soil. The roots were sampled before fertilizing and watering at 9 am, when the soil was dry and the seedlings were exposed to light 2 h. The materials collected from each plant were grounded to fine powder using liquid nitrogen and mixed, and send to BGI Life Tech Co., Ltd (Shenzhen, China) for Illumina sequencing.

RNA extraction, library preparation and sequencing
Total RNA was extracted using TRIzol 1 Reagent (Invitrogen) following the manufacturer's instructions. The integrity of the RNA samples was analyzed by agarose gel electrophoresis, and the purity and concentration were evaluated using an Agilent 2100 Bioanalyzer and Nano-Drop (Thermo Scientific, USA). All RNA samples had an RNA integrity number (RIN) higher than 8. Libraries were prepared using Illumina TruSeq RNA Sample Preparation Kit (Illumina, USA) according to the manufacturer's instructions. After checking quality and quantity, mRNA was isolated from total RNA using poly-(T) oligo-attached magnetic beads. The isolated mRNA was then sheared into short fragments by mixing with fragmentation buffer and used as templates for cDNA synthesis. cDNA was purified with a QIAquick PCR Purification Kit (QIAGEN), followed by end repair and acetylation of the 3' ends. The Illumina pair-end adapters were then connected with the cDNAs. Finally, the cDNA libraries were checked using a Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-time PCR system and sequenced using an Illumina HiSeqTM 2000 (Illumina Inc., San Diego, CA, USA). The sequence reads were submitted into the NCBI Short Read Archive (SRA) database and included in NCBI BioProject database with accession number PRJNA309919 (http://www.ncbi.nlm.nih.gov/bioproject/ 309919).

De novo assembly and gene annotation
After the filtering of low-quality reads and removing adaptors, clean reads were assembled using Trinity de novo transcriptome assembly [80]. The unigenes in each sample were processed for sequence splicing and redundancy removal to generate a non-redundant unigene dataset. The core eukaryotic gene-mapping approach (CEGMA) was used to assess the completeness of the de novo assembly. The unigenes were then searched against protein and nucleotide databases, including SwissProt, NCBI's non-redundant protein database (NR), NCBI's non-redundant nucleotide database (NT), Kyoto Encyclopedia of Genes and Genomes (KEGG) and COG (Cluster of Orthologous Groups) using BLASTX with an e-value cutoff of 1e-5. The unigenes were annotated by GO with Blast2GO sing the BLAST results for NR [81]. The unigenes were assigned with KEGG annotations using the BLAST results for KEGG [82]. Based on the results of protein database searches, the sequence direction and coding region were determined for unigenes with BLAST hits; ESTscan was used to predict coding regions for those without BLAST hits [83].

Unigene expression analysis
FPKM was calculated to assess unigene expression. NOISeq, a nonparametric approach, was used to determine the DEGs of samples with replicates, with a fold change2 and probability0.8 [34]. To characterize the function of DEGs, KEGG pathway enrichment analysis was performed using KOBAS 2.0 (q-value0.05) [83]. DEGs were also searched in PlantTFDB to determine DEGs related to TFs [35].

Validation of qRT-PCR
Fifteen DEGs were selected for qRT-PCR validation of gene expression. Gene-specific primers were designed with PRIMER 6.0 software (University of Plymouth); the primers are listed in S3 Table. Total RNA was extracted from SRs and FRs described as above. The total RNA were treated with DNaseI (Takara).With validation of the expression level of several reference genes using qRT-PCR and expression stability analysis the reference genes by Genorm (https:// genorm.cmgg.be) [84], the Actin-2 gene was used as an internal control (S3 Fig). The standard curve of each selected DEG was obtained by qRT-PCR with a series of cDNA dilutions. The 10 μl reaction mixture for qRT-PCR consisted of 5 μl of 2× SYBR Green Master Mix Reagent (Applied Biosystems), 0.2 μM of gene-specific primers and 50 ng of cDNA sample. The amplification reactions were performed with following program: 95°C for 10 min and 45 cycles of 95°C for 5 s and 60°C for 30 s. The relative expression levels of DEGs were calculated with the 2 −44Ct method. qRT-PCR was performed with 3 biological replicates and 3 technical replicates for each experiment.

Author Contributions
Conceived and designed the experiments: ZL LX ZW.
Performed the experiments: ZL ML MA YF.