Continuous inflorescence removal changes parameters, lignin and saponin a content and alters the transcriptome of root of Bupleurum chinense DC

We assessed the effect of the continuous inflorescence removal (CIR) on improving the yield and quality of Radix Bupleuri. The results showed that the taproot length, root head diameter, and lateral root numbers were significantly elevated by CIR. We also found an average 1.71-fold increase in root biomass, and the lower lignin and higher saponin a content were detected in roots of CIR-treated plants. Based on a comparative transcriptome analysis, 172, 243, 1974 and 3024 differentially expressed genes (DEGs) were respectively identified in the CIR-treated roots. Gene Ontology and KEGG pathway functional annotation analysis of the DEGs showed that multiple genes were involved sugars metabolism pathway, lignin and terpenoids biosynthesis pathway, and the plant hormones signal transduction pathway. In addition, 295 genes encoding transcription factors (TFs), including members of the ERF, MYB-related, bHLH, NAC, MYB, AUX/IAA and bZIP families etc, were identified as CIR responsive. first insensitive;

in rural areas of China. Nevertheless, the flowering period of cultivated B. chinense DC. could last as long as three or four months during a growing season. And previous study had found that because of the culturing condition, cultivated B. chinense DC. could produce high number of flowers and fruits, which seriously influenced the root yield [5]. The biomass of flower and fruits could reach to eight folds of the root biomass [6]. So it is meaningful and urgent to develop a method for high efficient cultivation of the roots of B. chinense DC.
In expection of decreasing the nutrient competition from reproductive growth and improve root growth, some investigators had examined the effects of inflorescence removal on plant root growth and chemical constituent content. Survey conducted by Kim had showed that the number of root, fresh root weight in peony were increased by two methods including removing flower buds or flowers, and the highest values were observed by removing flower buds [7]. Root biomass was increased when the inflorescence buds were partly removed in the plants of Panax ginseng C.A. Meyer, and total ginsenosides content of roots did not reduced [8]. Similar results of the increased biomass were also observed in ahipa [9], north American ginseng [10], potato [11], tomato [12], Helianthus tuberosus L. [13], and anchote [14] after inhibition of the reproductive growth. It was also showed that the operation of disbudding and picking flower increased the root biomass without effecting saikosaponin a and d contents of B. falcatum L. [15]. These finding suggested that inflorescence removal had a 4 positive effect on biomass production accumulation of roots or underground storage organs. However, flower bud thinning increased the root dry weight of "Misty", a variety of southern highbush blueberry, but not changed the root dry weight of "Santa Fe" [16]. On the other hand, although these studies has surveyed and proved the biomass increasing effect of inflorescence removal in some spieces, the related mechanism has not been investigated in any species ever.
In the present study, continuous inflorescence removal was conducted on plants of B. chinense DC.
through field experiments. The growth parameters, saponin, and lignin contents of roots were investigated. Moreover, the de novo transcriptome sequencing was performed to look for the possible molecular mechanisms underlying these changes of the root phenotypes.

CIR treatment changed the root growth parameters, lignification and saponin content
To investigate the influence of the CIR treatments on the market value of B. chinense DC. plants. The inflorescences in treated plants (defined as T) were removed as soon as inflorescence buds emerged since the day of the first inflorescence bud appeared. The treatment was continuously conducted for 75 days. Then roots were sampled for parameters measurement and chemical components analyses.
Root biomass was an important index of market value of cultivar of B. chinense DC.. The degree of lignification was also important for its market value since the dried roots of B. chinense DC. were cut into decoction pieces in traditional processing and serious lignification means higher mechanical strength and processing cost. Thus, we examined the effects of CIR on the root biomass and growth parameters, the total xylem cross-sectional area (CSA), total lignin contents, as well as the contents of saponin a, c, d which are the main medicinal ingredients of B. chinense DC..
The results showed that the biomass and growth parameters including root fresh weight, root dry weight, the taproot length, root head diameter and lateral root number of B. chinense DC. were all significantly higher in treated plants than in control plants. Especially for the fresh weight, dry weight, and lateral root number which showed 1.71-, 1.76-, and 0.87-fold increase, respectively ( Fig. 1a-e).
For the lignification investigation, although there was little difference of xylem CSA between treated and control plants, the lignin contents in roots of treated plants were found to be a 33.73 percent 5 decrease ( Fig. 1f and g). Thus the CIR treatment reduced the level of lignifications. The microstructures of root were also observed with light microscope. Like control plants (Additional file 1: Fig. S1a-d), intact anatomical structures in cross-sections of roots could be also observed in treated plants (Additional file 1: Fig. S1e-h). Furthermore, HPLC analysis showed that the content of saponin a was significantly higher in treated plants with a 32.01 percent increase, while the content of saponin c and d were not significantly changed (Fig. 1h-j).

Numbers of differentially expressed genes increased along with the duration of CIR
The differentially expressed genes (DEGs) were identified by pair-wise comparisons of these libraries, Obviously, the DEGs were obviously increased in the later stages.
In order to identify common and unique DEGs in response to CIR at different time points, a Venn diagram was plotted (Fig. 2b). Among these DEGs, none common DEG was identified in four groups.
Little common DEG genes were found for the first three time points. Four genes were shared between the "T1 vs FB" and "T2 vs IG" comparisons and 17 genes were shared between the "T2 vs IG" and "T3 vs R" comparisons. However, the "T3 vs R" and "T4 vs H" shared 738 DEGs.
To validate the accuracy of RNA-seq data, 16 DEGs were selected stochastically from transcriptome dataset and qRT-PCR analysis was carried out (Additional file 2: Table S1). Compared with the sequencing results, the qRT-PCR results showed the similar trends (Additional file 3: Fig. S2). Hence, the reliability of the expression level of RNA-seq data was confirmed.
Thus, the operation of CIR triggered or depressed the expression of various genes, and the DEG numbers of four groups displayed a markedly increasing tendency along with the duration of CIR operation.

GO and KEGG pathway enrichment analysis of DEGs
To further investigate the main biological functions of these DEGs from four groups, the DEGs were subjected to the GO and KEGG database. In general, 4813 DEGs from four groups were GO-annotated into 44 terms. 172 DEGs in "T1 vs FB" were GO-annotated into 31 terms. In "T2 vs IG", 243 DEGs were GO-annotated into 30 terms. In "T3 vs R", 1974 DEGs were GO-annotated into 43 terms, and 3024 DEGs in "T4 vs H" were GO-annotated into 41 terms ( Fig. 3 and Additional file 4: Table S2).
Although the DEG numbers were diverse, "cellular process", "metabolic process", and "singleorganism process" were the most enriched terms in biological process (BP), and "catalytic activity" and "binding" were the most enriched terms in molecular function (MF) for all the four time points.
There were five, four, six and five KEGG categories based on pathway hierarchy1 (Additional file 5: Table S3). Among the pathway subgroups based on pathway hierarchy 2, the pathway of Carbohydrate metabolism contained the most genes, followed by Lipid metabolism (Fig. 4 and Additional file 5: Table S3). All of the KEGG-annotated DEGs were assigned to 23, 31, 86 and 97 KEGG pathways, respectively (Additional file 5: Table S3). For DEGs of "T1 vs FB" and "T2 vs IG", "sesquiterpenoid and triterpenoid biosynthesis" and "Protein processing in endoplasmic reticulum" was the top enriched pathway, respectively. The top four enriched pathways in both "T3 vs R" and "T4 vs H" were "starch and sucrose metabolism", "pentose and glucuronate interconversions", "plant hormone signal transduction" and "phenylpropanoid biosynthesis" (Additional file 5: Table S3) .

Major DEGs profiling of sugar metabolic pathways
Based on the KEGG pathway enrichment analyses, DEGs involved in the carbohydrate metabolism were the most abundant, so we analyzed the "starch and sucrose metabolism", "galactose metabolism" and "pentose and glucuronate interconversions" pathways. 38 and 48 DEGs from "T3 vs R" and "T4 vs H" were enriched in "starch and sucrose metabolism" pathway, respectively. Among them, 36 DEGs were up-regulated and two were down-regulated in "T3 vs R", while 38 DEGs were up-regulated and 10 were down-regulated in "T4 vs H" (Fig. 5a and b and Additional file 6: Table S4). Four and five DEGs from "T3 vs R" and "T4 vs H" were enriched in "galactose metabolism" pathway, respectively (Additional file 6: Table S4 and Additional file 7: Fig.   S3a and c). Both "T3 vs R" and "T4 vs H" had 31 DEGs enriched in "pentose and glucuronate interconversions" pathway (Additional file 6: Table S4 and Additional file 7: Fig. S3b and d). Most of 7 these DEGs were up-regulated. In "T3 vs R", 29 genes were up-regulated, while in "T4 vs H", the number was 28. The annotation showed that some DEGs including xylosidase gene, xylan synthase gene, UDP-glucuronate decarboxylase gene (UXS), α-galactosidase gene, β-fructofuranosidase gene, endoglucanase gene and glucosidase gene were related to the xylose, galactose and glucose biosynthesis and metabolism (Additional file 6: Table S4 Table S4).

Expression of DEGs involved in phenylpropanoid and terpenoid biosynthesis pathways
Lignin was generated through phenylpropanoid biosynthesis pathway which had been well elucidated [17]. Saponin was biosynthesized through terpenoid biosynthesis pathway [18]. In our study, the DEGs from phenylpropanoid and terpenoid biosynthesis pathways were detected. Three, 19 genes and 34 genes from "T2 vs IG", "T3 vs R", and "T4 vs H" were enriched in phenylpropanoid biosynthesis pathway, respectively. Among the DEGs, some upstream genes were up-regulated. For instance, one phenylalanine ammonia lyase gene (PAL) was up-regulated in "T4 vs H", and three hydroxycinnamoyl-CoA:shikimate/quinate hydroxycinnamoyltransferase genes (HCT) were upregulated in "T4 vs H". One and six beta-glucosidase genes were up-regulated by CIR in "T2 vs IG" and "T3 vs R", respectively. One caffeic acid O-methyltransferase gene (COMT; Cluster-10318.50166) was up-regulated in "T3 vs R", and another gene in this family (Cluster-10318.33017) was upregulated in "T4 vs H". Interestingly, we observed the serious reduction of the expression level of seven genes encoding cinnamyl alcohol dehydrogenase (CAD) in "T3 vs R" or "T4 vs H" group. CAD is 8 a key enzyme in monolignin biosynthesis, and down-regulation of CAD gene generally leads to reduction of lignin content [19]. The last step in lignin biosynthesis pathway is catalyzed by peroxidase (POD) which transforming cinnamyl alcohols to lignin. A number of POD genes were identified in this study. Both up-regulated genes and down-regulated genes were exhibited in this gene family. (Fig. 6 and Additional file 8: Table S5).
Likewise, CIR modulated the expression of most terpenoid biosynthetic genes in four groups, and most of these genes were up-regulated in the last three groups. Particularly, these DEGs including  Table S6).

9
The expression of genes involved in hormone signal transduction pathways was altered following the CIR treatment. According to the KEGG analysis of DEGs, 30 and 46 DEGs in this pathway were identified from "T3 vs R" and "T4 vs H", respectively (Fig. 8). Among different hormone signal transduction pathways, auxin signal transduction pathway contained the largest number of identified DEGs for both the "T3 vs R" and "T4 vs H" groups, followed by the cytokinine and brassinosteroid signal transduction pathways.
Auxin, inducing three kinds of early auxin response genes: AUXIN/INDOLE-3-ACETIC ACID (AUX/IAA), GRETCHEN HAGEN3 (GH3), and SMALL AUXIN UP RNA (SAUR) which had been identified in diverse plant species [20], was supposed to be the primary hormones that regulate the lateral roots development and elongation [21]. We observed significant changes of 10 and 11 genes involved in auxin signal transduction pathways of "T3 vs R" and "T4 vs H", respectively, including one gene encoding transport inhibitor response protein (TIR), twelve auxin-responsive protein genes (nine AUX/IAA genes and three SAUR genes), and five auxin response factor genes (ARF). TIR functioned as auxin receptor [22]. Most of AUX/IAA genes were up-regulated, while all ARF genes were downregulated in "T3 vs R" and "T4 vs H". In addition, the three genes encoding SAUR included two upregulated genes and one down-regulated gene, and one gene (Cluster-10318.15908) had a 3.39-fold change in "T4 vs H" (Additional file 11: Table S8).
Among the cytokinine signal transduction-related genes, 100% of DEGs (6/6) were up-regulated in "T3 vs R", and 66.7% of DEGs (6/9) were up-regulated in "T4 vs H". Three and four genes encoding histidine kinases (HKs), the cytokinine-binding receptors, were up-regulated in "T3 vs R" and "T4 vs H", respectively. And one HK gene was down-regulated in "T4 vs H". A multistep HK→histidinecontaining phosphotransfer protein (HP)→response regulator (RR) phosphorelay pathway was involved in response to cytokinine [23]. One gene encoding HP functioned as phosphor-histidine intermediates between the HK and RR was up-regulated in "T4 vs H". Three RR genes were up-regulated in "T3 vs R", and while one RR gene was up-regulated and two RR genes were down-regulated in "T4 vs H" (Additional file 11: Table S8).
In brassinosteroid (BR) signal transduction pathway, 80% of DEGs (4/5) were up-regulated in "T3 vs R", and 77.8% of DEGs (7/9) were up-regulated in "T4 vs H". Brassinosteroid insensitive (BRI) functioned as the brassinosteroid receptor and it is a member of the leucine-rich repeat receptor-like kinase family. One gene encoding BRI kinase inhibitor which is the inhibitor of brassinosteroid receptor (BKI) was down-regulated in "T3 vs R" and "T4 vs H". Moreover, two brassinosteroidsignaling kinase genes (BSK, Cluster-10318.45934 and Cluster-10318.84972), the important positive regulator of BR signal, were up-regulated in "T3 vs R" and "T4 vs H". In general, BR triggered active cell division associated with increased transcripts of cell cycle-related genes, especially that of cyclin D3 genes [24]. Notably, seven cyclin D3 genes were identified and all genes were up-regulated by CIR treatment (Additional file 11: Table S8).

Discussion
In previous works, the root biomass or length could be up-regulated after removal of the flower buds, floral or inflorescences in ahipa [9], north American ginseng [10], potato [11], tomato [12], Helianthus tuberosus L. [13], anchote [14], B. falcatum L. [15], Platycodon grandiflorum A. DC. [25], Duchesnea indica [26] and cotton [27]. However, secondary metabolites in roots were not changed by the treatment in the majority of these researches. Our data showed that the CIR treatment could not only increase root biomass and growth parameters, but also elevate the content of saponin. This might be caused by the increase of treat time and degree in our study. Furthermore, we found that CIR reduced the lignin content significantly, which could provide convenience for root processing. These finding indicated that CIR operation had great potential to improve the yield and quality of Radix Bupleuri.
Furthermore, this approach may also be conducted in the cultivation of other plants which provide roots as traditional medicine or food to improve their yield, quality, and market value.
Because of the importance of CIR treatment in cultivated B. chinense DC. and its possible importance in other medicinal herbal with root as medicinal part, which account for a large proportion of traditional Chinese medicine, and root vegetables, we monitored the overall gene expression level in roots of B. chinense DC. during CIR treatment by RNA-seq. We found that the DEGs number gradually increased with the duration of CIR, and dramatically increased in latter two groups, it suggested that the changes in gene expression and possible root phenotype differences were enlarged with the treatment time lasting.
Our results found that almost all DEGs of three carbohydrate metabolism pathways ("starch and sucrose metabolism", "galactose metabolism" and "pentose and glucuronate interconversions") were up-regulated in "T3 vs R" and "T4 vs H", and these DEGs were related to xylose, glucose, and galactose metabolism. For example, The UXS genes which were up-regulated here were positive regulators of xylose, glucose, and galactose in tobacco or Arabidopsis by the RNAi transgenic experiments [28,29]. These genes' up-regulation indicated the activation of the production of xylose, galactose and glucose related metabolites. Since these sugars are important in the biosynthesis of saponin [2], these DEGs might play important roles in the increase of the saponin content after CIR treatment. Until now, no study has presented direct effects of these DEGs on the saponin content.
Thus, digging of the functions of these genes in saponin production will be an interesting direction for future study. On the other hand, we found a lot of DEGs were annotated as endoglucanase genes and PL genes. Endoglucanase (EC 3.2.1.4) was an enzyme required for biosynthesis of cellulose, and it had been proved that it affected the root elongation. The endo-1,4-beta-glucanase mutant of rice exhibited defects in both root cell elongation and division, and had lower cellulose content, shorter root and lateral root than wild type [30]. Pectin was a major component of primary plant cell wall. PL (EC4.2.2.2), a cell wall loosening/hydrolytic enzyme [31], catalyzed the breakdown of pectin of cell walls in plant organs. Main root length of PL mutant of rice was clearly reduced and the number of lateral roots significantly decreased, along with the delayed cell cycle progression and less cell number [32]. It is possible that the up-regulation of these genes jointly promoted cell division and elongation, thus the root parameters changed.
Since we found the lignin content was reduced after CIR treatment, we further investigated the expression of genes in phenylpropanoid pathway and found the significantly reduction of CAD genes' expression. CAD is a key enzyme in monolignin biosynthesis, down-regulation of CAD genes generally leads to reduction of lignin content [19]. Peroxidases were encoded by multigene families and involved in wide range of physiological processes such as lignification, auxins metabolism, wound healing and reactive oxygen species (ROS) in plants [33]. A higher level of lignin was obtained in the transgenic tobacco with overexpressed peroxidase gene than in the wild-type plants [34]. Here, both up-regulated and down-regulated peroxidase genes were identified in "T2 vs FR", "T3 vs R" and "T4 vs H". These down-regulated genes might participate in the lignin biosynthesis of B. chinense DC.
root. Thus, we supposed that these down-regulated CAD genes combined with peroxidase genes might be reasons for less lignin content.
The majority of transcripts associated with terpenoid biosynthesis were up-regulated in this study.
Saikosaponins (triterpenoid saponins) were the main effective components of B. chinense DC., and formed through terpenoid backbone biosynthetic pathway including MVA and MEP pathways [18]. Several TFs (such as AP2/ERF, MYB, and GRAS) required for root development were found in our DEG data. Among up-regulated AP2/ERF genes in "T3 vs R", one was homologous to CRF2 gene (Cluster-10318.91577), a positive regulator that increased lateral root densities when it was overexpressed [40]. The plant with overexpressed MYB59 gene suppressed root growth, and loss of function MYB59 mutant had longer primary root than wild type plant [41]. Here, all differently expressed MYB59 genes homologies were repressed after CIR in "T3 vs R" and "T4 vs H".
Lateral root formation relied on the generation of lateral root primordia (LRP). One up-regulated TF gene was annotated as MYB36 gene, and MYB 36 has been shown as a key regulator in the transition from proliferation to differentiation in the primary root [42], and was also necessary for LRP emergence through activating expression of peroxidase gene [43]. Moreover, three SCARECROW-like 28 genes belonging to GRAS genes were up-regulated with ~3 fold change in "T4 vs H". SCARECROW was proved that could directly activate the expression of MYB36 genes [42]. These data indicated root formation and development after CIR treatment was directly or indirectly regulated by these TFs. The detailed mechanism remained to be elucidated.
In this study, we found the three SAUR genes and two of them were up-regulated after CIR treatment.
Overexpression of SAUR 41 gene in Arabidopsis led to auxin-related phenotypes including increased vegetative biomass and lateral root development [44]. Moreover, one TIR1 gene was also upregulated. The number of lateral roots was reduced in TIR1 mutant of Arabidopsis seedlings [45]. So it suggested that auxin might play an important role in the root phenotype alteration after CIR. BR also plays a crucial role in regulating many of the aspects of plant growth and development, including cell division, expansion and root growth [46]. Intriguingly, all identified DEGs that encoded brassinosteroid-signaling kinase were up-regulated. BKI1 is negative regulator of BRI1 function [47].
Furthermore, the expression of a BKI gene was down-regulated in "T3 vs R" and "T4 vs H". Based on these results, we proposed that the BR signal pathway might be activated in treated plants of B.
chinense DC. and lead to the alteration of root growth. The cytokinine was mainly responsible for regulating cell division and differentiation [48]. We detected the up-regulated DEGs related to HK and RR. HK 4 was primary receptor kinases that could transduce cytokinine signals across the plasma membrane in Arabidopsis thaliana [49]. And in Arabidopsis thaliana, ARR5 and ARR17 were immediately induced by cytokinine at the level of transcription [50]. The primary root length of overexpression of ARR5 and ARR17 in Arabidopsis plants was significantly longer than that of wild type [51]. Moreover, we found DEGs that encoding cyclin D3 were all up-regulated by CIR. The transcription of cyclin D3 as a regulatory protein of the cell cycle was positively enhanced by BRs in previous report [52], and the expression of cyclin D3 transcript was also regulated by cytokinine [53].
This demonstrated that BR and cytokinine might influence the growth of B. chinense DC. root through 14 cyclin D3.

Conclusion
In the present study, we found that the biomass, growth parameters and saponin a content were elevated, but the lignin content was decreased by CIR.. These results demonstrated that CIR treatment could improve the yield, quality and market value of B. chinense DC. Moreover, this study is the first report of transcriptome data of the roots in plants under CIR condition. The present study highlighted several key pathways and multiple genes that might play important roles in CIR response, these genes were involved in sugars metabolism, lignin and terpenoids biosynthesis. 295 genes encoding TFs, including members of the ERF, MYB-related, bHLH, NAC, MYB, AUX/IAA and bZIP families etc, were identified as CIR responsive. Our results will facilitate uncovering the molecular mechanism of root when responding the CIR and provide a valuable method to promote biomass and the content of chemical components in the root herbs. Plants were established at a density of 9 plants/m 2 . All cultivation measures were applied according to the regional recommendation. When all plants grow up adequately, they were randomly divided into two groups, the treated group (TG) and the control group (CG) and tagged. The inflorescences were removed at one or two days intervals using scissor as soon as inflorescence buds emerged at the earliest date that was July 29, 2018 (Additional file 12: Fig. S4a), and the operations were persistent till October 11, 2018, in TG. Roots were sampled at four time points during the 75 days' treatment: one in the full bloom stage (FB stage, Additional file 12: Fig. S4b), one in fruit immaturegreen stage (IG stage, Additional file 12: Fig. S4c), one in fruit ripeness stage (R stage, Additional file 15 12: Fig. S4d), and one in root harvest stage (H stage, Additional file 12: Fig. S4e). Finally, a total of 24 samples including four group pairs of treated and control samples ( T1 vs FB, T2 vs IG, T3 vs R, and   T4 vs H ) were collected. The fresh roots for RNA sequencing were separated and washed carefully with double distilled water, then were tightly packed in tinfoil and immediately immerged in liquid nitrogen and stored at -80°C prior to use. At the last sampling date, all plants were collected, and divided into several subgroups for determining root morphology and the contents of chemical components.

Root Growth Measurements, Anatomy and Chemical Analyses
At H stage, fifteen roots in TG and CG collected in October 11, 2018 were assigned to the following experiments. The maximum lengths, the fresh weights, the head diameters of the intact root system freshly excavated were measured by ruler, analytical balance and vernier caliper, respectively.
Number of lateral roots (diameter≥2mm) was counted. Fresh roots were placed in a freeze dryer until complete drying and dry weights were determined.
Three roots from TG and CG were randomly selected and root segments (approximately 0.5cm) were hand-made with a razor blade and stored in 50% FAA until further processing. The roots were The significance between treated and control group was determined using Student's t-test.

RNA Extraction and Deep Sequencing
The total RNA from 24 root samples collected in four dates as mentioned above were extracted using All sequencing libraries were generated using NEB Next Ultra RNA Library Prep Kit for Illuminu (NEB, USA) following manufacturer's recommendations. The literary quality was assessed on the Agilent Bioanalyzer 2100 system. In brief, mRNA was purified from total RNA and fragmentation was carried out. First strand and second strand cDNA synthesis were in turn performed using random hexamer primers and reverse transcriptase. The library preparations were sequenced on an Illumina Hiseq 2100 platform (Novogene, Beijing, China). Raw data of FASTQ format were firstly processed through in-house Perl scripts. Reads containing adapter, reads containing ploy-N and low quality reads from raw data were removed and clean reads were obtained. In order to assess the sequencing quality, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality. Transcriptome assembly was accomplished using Trinity v2.4.0 with min_kmer_cov set to 2 by default and all other parameters set default [56].

Transcript Annotation and Gene Expression Analysis
Gene function was annotated based on the following seven databases: To quantify expressions levels of the genes, the read counts were calculated by using the fragments per kilobase of transcript per million fragments mapped method (FPKM) [57].

Identification of Differentially Expressed Genes
Differential expression analysis of different groups was performed using the DESeq2 R package (version 1.6.3). The adjusted P-value (padj) <0.05 and an absolute log2 fold change>1 was set as the threshold for significantly differential expression. Venn diagram analysis of DEGs between different groups was performed with an online tool (https://magic.novogene.com/public/customer/login).
Heatmaps of different gene expressions level were exported using the software of HemI (version 1.0, http://hemi.biocuckoo.org/). Gene Ontology (GO) enrichment analysis of DEGs was implemented by the GOseq R packages based Wallenius non-central hyper-geometric distribution, which can adjust for gene length bias in DEGs [58]. We used KOBAS software to test the statistical enrichment of differential expression genes in KEGG pathways [59].

Availability of data and materials
The data sets generated or analyzed during this study are included in this published article and its additional files. All the transcriptome data used and analysed during the current study are available from the corresponding author on reasonable request.

Ethics approval and consent to participate
Not applicable.

Consent for publication
Not applicable.

Competing interests
The authors declare no conflict of interest.  DEGs assigned to "starch and sucrose metabolism" pathways in "T3 vs R" (a) and "T4 vs H" (b) under CIR condition. Red box up-regulated genes; cyan box down-regulated genes; yellow box up-regulated and down-regulated genes.