Effect of different pruning operations on the plant growth, phytohormones and transcriptome profiles of the following spring tea shoots

Tea, produced from the tender shoots of Camellia sinensis (L.) O. Kuntze, is a beverage consumed worldwide. Tea is associated with a variety of health benefits, such as anti-inflammatory^[1], anti-tumorigenic, and cardioprotective effects^{[2, 3]}. The quality of fresh tea leaves is important to the quality of final tea products, which is affected by several factors, such as harvest season^[4], cultivation practice^[5], climatic conditions^[6], and plucking standard. Pruning is one of the conventional agronomic operations for tea cultivation^{[7, 8]}. Pruning treatments can regulate the growth of branches and young shoots^[9], shape the tea trees, and affect the production yield and the quality of fresh tea leaves^{[9, 10]}. The accumulation of secondary metabolites in tea leaves was also affected by pruning operations^{[7, 11]}.

Different conventional pruning operations are employed in tea cultivation practice for different purposes. There are light pruning, deep pruning, heavy pruning and collar pruning with increasing pruning depth. Specifically, heavy pruning refers to cutting the main stems back to 40−45 cm from the ground. Deep pruning means cutting the main stems by 10−15 cm from the top of the tea trees. Light pruning means cutting the main stems by 3−5 cm from the top of the tea trees. Heavy pruning and collar pruning with greater pruning depth usually exert more effective regenerative effect on tea branches for renewing or removal of diseased branches than light pruning and deep pruning^[12], while the production yield is greatly reduced in the following year after heavy pruning and collar pruning^{[13, 14]}. Mohale et al. reported that cutting the main stems by 15 cm from the top caused greater accumulation of metabolites in the tea leaves after 90 days, compared with other treatments like heavy pruning and collar pruning^[15]. Jiang et al. reported that selective pruning in summer improved the production yield and quality of spring tea compared with deep pruning^[16]. Pruning at a height of 45 cm from the ground is helpful to form the plucking surface of tea plants. Moreover, pruning time is also an important factor in the growth of tea plants, especially for the production yield and quality of the following spring tea^[9,17]. It was reported that light pruning in late spring could improve the production yield and quality of the following spring tea shoots, compared with light pruning in winter^[18]. Summer pruning (cutting the main stems by 30 cm from the top) is more effective to promote the growth of lateral branches as well as the production yield and quality of spring tea, compared with autumn pruning^[9]. It was reported that in India, the higher yield was obtained as pruned in December, compared with pruning in November, which could be associated with the high accumulation of root starch in December^[19]. Deep pruning after spring tea is suitable for producing green tea, while deep pruning after summer tea is suitable for producing black tea^[20]. Besides, the flavor index value for tea increased the year after pruning, which was consistent with the results of sensory evaluation^[10]. In a 5-year pruning cycle, the production yield reached a peak the 3^rd year after pruning, and then decreased^[19]. Due to the effects of pruning depth and time, the combination of pruning depth and time may have comprehensive effects on the growth and metabolism of young shoots. However, most studies focused on the production yield and chemical composition of fresh tea leaves under different pruning operations, and the development stage of young tea shoots crucial to production yield and chemical composition of fresh tea leaves was not fully taken into account. There is still a lack of systematic studies on the impacts of pruning depth and pruning time on the growth and metabolism of the following spring tea shoots.

The growth rate of spring tea shoots is closely related with the economic benefits of tea production, since the early spring tea of premium quality can be sold at a higher price^[21]. Here, the effects of different pruning operations on the growth of tea plants were investigated, including different pruning time [middle April (4.15), early May (5.8), and late May (5.28)] and different pruning depth [light (L), deep (D) and heavy (H) degree]. The profiles of phytohormones, catechins and transcriptomes in the following spring buds after different pruning treatments were studied. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was carried out to reveal the significant biological changes of the following spring buds after different pruning treatments, and the genes involved in the enriched pathways of plant-pathogen interaction, plant hormone signal transduction and circadian rhythm were further interpreted. This study provides scientific guidance for tea plant pruning operations, and preliminarily reveals the regulatory mechanism of pruning on the growth of spring shoots from the aspects of phytohormone and transcriptional levels.

The pictures of spring shoots and buds in the following year after different pruning treatments are shown in Fig. 1. Obviously, heavily pruned tea plants had thick and sparse branches, which were quite different from the control (CK), while the branch morphology of tea plants after light and deep pruning was in between. The buds of HP 4.15, HP 5.8 and HP 5.28 were obviously bigger and stronger than those of CK, DP 5.8, LP 5.8 and HP 5.8 + LP 9.30 (Fig. 1). Figure 2 shows the development stages of the following spring tea shoots. The tea shoots of HP 4.15 and HP 5.8 grew faster than that of CK, while the development of the tea shoots of HP 5.8 + LP 9.30 and HP 5.28 were comparable with CK. For example, on March 15^th, the young shoots of HP 4.15 and HP 5.8 were at the stage of two leaves and one bud, while the young shoots of CK and other pruning treatments were at the stage of one leaf and one bud, suggesting that the tea plants that were heavily pruned in the middle April and early May had early sprouting of young shoots in the following spring, compared with CK and other pruning treatments.

Figure 1. 

The effects of different pruning treatments on the spring shoots and buds in the following year.

Figure 2. 

The developmental stage of young tea shoots during the following spring after different pruning treatments. The early development of young tea shoots can be divided into five characteristic stages: sprouting, first fish leaf expansion, second fish leaf expansion, one leaf and one bud, as well as two leaves and one bud.

The box plots of stem diameter, shoot length and weight of 100 buds (two leaves and one bud) for the following spring tea shoots with different pruning treatments are shown in Fig. 3. At the stage of two leaves and one bud, the median of the stem diameter for the pruned tea plants ranged around 1.8 mm, which was higher than that of CK (1.7 mm). Moreover, the average value of the stem diameter for the pruned tea plants ranged from 1.7 to 1.8 mm, compared with 1.7 mm of CK. For shoot length, the median and average value of CK were 87.4 mm and 84.7 mm, which were much higher than that of the pruned tea plants. By contrast, the median and the average value of the weight of 100 buds were 14.0 g and 14.0 g for CK, which were much lower than those of the pruned tea plants in the range of 15.5−20.5 g and 16.7−19.9 g. The relatively lower weight of 100 buds but longer shoot length for CK means the buds of CK were much thinner than those of pruned tea plants.

Figure 3. 

The box plots of (a) stem diameter, (b) shoot length and (c) weight of 100 buds of young tea shoots in the following year after different pruning treatments. Stem diameter was measured in the middle of the third internode from the apical bud; shoot length refers to the length from the base of stem to the top of apical bud; weight of 100 buds was weighted at the stage of one bud and two leaves. All the parameters were measured at the stage of two leaves and one bud on March 22^nd, 2021.

Effects of different pruning operations on the levels of phytohormones in the following spring buds

Phytohormones importantly regulate the growth of young tea shoots during the sprouting period. Table 1 shows the levels of endogenous phytohormones at the bud stage after different pruning treatments. The pruning-treated buds generally contained higher phytohormone levels than CK. Moreover, HP 5.8 contained the highest levels of auxin (IAA), gibberellin A1 (GA1), trans-zeatin (tZ), gibberellin A3 (GA3) and trans-zeatin riboside (tZR), subsequently followed by HP 5.8 + LP 9.30 and DP 5.8, while CK and HP 5.28 contained low levels of these phytohormones. HP 4.15 had the highest level of jasmonic acid (JA) at 6,399.88 pg/g dry weight (DW), while HP 5.8 + LP 9.30 contained the highest level of melatonin (516.89 pg/g DW). Thus, pruning significantly affected the phytohormone levels of tea buds from the aspects of pruning time and pruning depth. For heavy pruning, the content of JA in buds decreased from 6,399.88 pg/g to 2,191.14 pg/g as the pruning time was postponed from 15^th April to 28^th May, while abscisic acid (ABA) showed the opposite change trend. For different pruning depth on the same day, the levels of GA1, GA3, tZ and tZR increased as the pruning depth increased, and HP 5.8 contained higher content of IAA, compared with DP 5.8 and LP 5.8.

Effects of different pruning operations on the content of catechins in the following spring buds

Table 2 shows the content of catechins in the buds of the following spring after different pruning treatments. HP 4.15 contained the highest level of total catechins (TC) at 58.25 mg/g DW, subsequently followed by HP 5.8 (56.24 mg/g DW), while no significant difference of TC was observed among CK, DP 5.8, LP 5.8, HP 5.8 + LP 9.30, and HP 5.28. As the most abundant catechin compound, EGCG showed the same change trend as TC. The different level of TC in buds might be related with the different growth stages of buds. HP 4.15 and HP 5.8 had an earlier development compared with CK and other pruning-treated tea plants.

Transcriptome profiles of differently pruned tea samples and KEGG enrichment analysis

Supplemental Table S1 shows the information of the RNA-Seq data. There were 41.6–43.4 million, 38.1–46.5 million, 38.1–55.0 million, 47.2–59.0 million, 51.1–59.1 million, 50.8–61.3 million, 48.4–51.2 million RNA-seq clean reads respectively obtained for CK, HP 4.18, HP 5.8, DP 5.8, LP 5.8, HP 5.8 + LP 5.8, and HP 5.28 (Supplemental Table S1). The mapping rates of reads to the reference genome were above 90% (analyzed by HISAT2). The PCA score plot of RNA-seq data is shown in Fig. 4a, and the first two principal components (PC) accounted for 64.9% of the total variance (PC1 = 50.8%, PC2 = 14.1%). Only HP 4.15 was close to CK, both of which were clearly discriminated from other pruning-treated samples at the direction of PC1. Moreover, the samples heavily pruned in May, including HP 5.8, HP 5.8 + LP 9.30, and HP 5.28, were mainly located at the positive direction of PC2, while DP 5.8 and LP 5.8 were located at the negative direction of PC2. The results of gene expression obtained by RNA-seq were validated by RT-qPCR (quantitative real-time polymerase chain reaction analysis), with R² being 0.8313 (Supplemental Fig. S1). Thus, transcriptome data were highly correlated with the RT-qPCR results, which could well represent the transcript abundances. Figure 4b shows the number of up-regulated and down-regulated DEGs in the different pruning-treated tea samples compared in pairs. The pairs of CK vs. HP 5.8 had the highest number of DEGs (114 up-regulated DEGs and 533 down-regulated DEGs), subsequently followed by the pairs of CK vs. HP 5.28 (99 up-regulated and 487 down-regulated DEGs) and the pairs of CK vs. HP 5.8 + LP 9.30 (120 up-regulated and 448 down-regulated DEGs). Less than 20 DEGs was observed in the pairs of HP 5.8 vs. DP 5.8, HP 5.8 vs. HP 5.8 + LP 9.30, HP 5.8 vs. HP 5.28, LP 5.8 vs. DP 5.8, LP 5.8 vs. HP 5.28, DP 5.8 vs. HP 5.28, and HP 5.28 vs. HP 5.8+LP 9.30, suggesting that these comparison pairs were highly similar in their transcriptional profiles.

Figure 4. 

The transcriptome profiles of the buds in the spring of following year after different pruning treatments. (a) PCA analysis, (b) DEG numbers, (c) significantly enriched KEGG pathways of DEGs. The number of replicates is three.

KEGG enrichment was employed to analyze the enriched pathways of DEGs. The top 10 pathways are shown in Fig. 4c. The comparison pairs of CK vs. pruning-treated samples as well as HP 4.15 vs. other pruning treated samples had more enriched pathways than the comparison pairs of pruning-treated samples in May, including RNA transport, plant hormone signal transduction, plant–pathogen interaction, photosynthesis–antenna proteins, circadian rhythm–plant. No significantly enriched pathway was observed in the comparison pairs containing the tea samples pruned in May. Specifically, the comparison pair of CK vs. HP 4.15 had the highest amount of significantly enriched pathways, including plant hormone signal transduction, MAPK signaling pathway, RNA transport and plant-pathogen interaction, subsequently followed by the pairs of CK vs. HP 5.8 + LP 9.30, HP 4.15 vs. HP 5.8 as well as HP 4.15 vs. HP 5.8 + LP 9.30 that were commonly enriched in the pathways of circadian rhythm–plant, photosynthesis–antenna proteins, and plant hormone signal transduction. Thus, pruning could differentially regulate the plant hormone signal transduction in tea leaves, which might affect the development of young tea shoots.

The expressions of DEGs in the pathways of plant–pathogen interaction, plant hormone signal transduction and circadian rhythm

Based on the result of KEGG enrichment analysis, the heatmaps of DEGs in the pathways of plant–pathogen interaction, plant hormone signal transduction and circadian rhythm are shown in Fig. 5. There were 59 DEGs obtained in the pathway of plant–pathogen interaction (Fig. 5a). Different samples were clustered into two major groups, including CK and pruning-treated samples based on the expressions of these 59 DEGs in the pathway of plant–pathogen interaction, and then the pruning-treated samples were sub-grouped into HP 4.15 and the samples pruned in May. Moreover, HP 5.8 and HP 5.8 + LP 9.30 had a high similarity, while DP 5.8 and LP 5.8 had a high similarity, which was close to HP 5.28. Differently, HP 5.8 + LP 9.30 and HP 5.28 had a high similarity of gene expressions in the pathway of plant hormone signal transduction, which were close to HP 5.8 (Fig. 5b). In Fig. 5c, CK and HP 4.15 had a high similarity of gene expressions in the pathway of circadian rhythm, which were differentiated from two subgroups. One subgroup consisted of HP 5.8 and HP 5.8 + LP 9.30, the other subgroup consisted of HP 5.28, DP 5.8 and LP 5.8. Thus, the circadian rhythm of buds was more affected by pruning time as heavily pruned, while plant–pathogen interaction and plant hormone signal transduction were more affected by pruning depth.

Figure 5. 

The heatmaps of DEGs in the fresh tea leaves under different pruning treatments. (a) Plant–pathogen interaction, (b) plant hormone signal transduction, and (c) circadian rhythm.

Pruning operation affects the production yield and the quality of fresh tea leaves^{[9, 18]}. Different results have been reported about the effect of pruning operation on tea plants^{[7−20, 22]}. Light pruning after spring tea increased the production of the following spring tea^[18]. Pruning in summer increased the number of branches per m², which was propitious to the production yield and quality of the following spring tea^[9]. However, the effect of pruning on the growth of young tea shoots is related with harvest time and control sample used. In the present study, heavy pruning in mid April and early May accelerated the development of the following spring tea shoots and elevated the weight of 100 buds at the stage of two leaves and one bud. The pruned tea plants generally had longer stem diameter and shorter shoot length, compared with CK, which was consistent with the previous study that the growth of new stems were thicker and shorter after heavy pruning treatment^[9]. The increased weight of 100 buds but reduced length of young shoots in the pruned tea plants indicated the weight ratio of leaves in the young shoots were much higher than that of the control sample. This possibly contributes to the generally considered elevated quality of fresh tea leaves in the following spring due to heavy pruning.

Various phytohormones are involved in the development of tea plants^[23]. Our study shows the levels of growth promoting phytohormones, such as IAA, GA1, GA3, and tZ, were significantly elevated in the spring buds of tea plants heavily pruned in early May, which might be a plausible explanation for its earlier development and better growth status. Moreover, the composition of catechins in tea leaves can be affected by pruning treatment^[22]. In the present study, the percentage of EGCG and ECG in TC ranged from 89.1% to 90.1% for the tea plants heavily pruned before the 8^th of May, which was higher than 88.8% of CK. This is in agreement with a previous study that the ratio of galloylated catechins were elevated in the leaves of pruned tea plants, compared with unpruned tea plants^[11]. Compared with pruned tea plants, more glucose, sucrose, catechin, EC, EGC, ECG, GC and epicatechin 3-O-(3-O-methyl)-gallate were accumulated in the leaves of unpruned tea plants^[7]. Besides, only HP 4.15 and HP 5.8 contained significantly higher contents of TC than CK, DP 5.8, LP 5.8, HP 5.8 + LP 9.30, and HP 5.28, while no significant difference of TC content was observed among these samples. The higher content of TC in HP 4.15 and HP 5.8 could be related to the earlier development of buds, since the composition of catechins is related with the growth stage of tea leaves^[24].

Pruning is a typical mechanical damage to tea plants, which arouses or regulates various physiological activities. In the present study, the pathways of plant–pathogen interaction, plant hormone signal transduction and circadian rhythm were significantly enriched in the following spring buds of pruned tea plants, based on the results of KEGG. The plant–pathogen interaction was involved with plant hormone signal transduction^[25]. The significantly enriched pathway of plant hormone signal transduction in the spring buds with heavy pruning treatment was in agreement with the regulatory effect of pruning operation on the levels of phytohormones in tea leaves. Compared with CK, the expression levels of most genes related to plant-pathogen interaction and plant hormone signal transduction were generally down-regulated in the pruning-treated samples. WRKYs and MYBs are important transcription factors playing crucial roles in plant responses to various abiotic stresses^[26,27]. The expressions of WRKY33 (TEA002467.1, TEA028505.1) and MYB44 (TEA018997.1) were remarkably down-regulated in the pruning-treated samples, compared with CK. In Arabidopsis, WRKY33 is a key transcriptional regulator of hormonal and metabolic responses toward fungal infection^[28]. MYB44 is a negative regulator of ABA, which appears to have a homeostatic function to maintain growth processes under stress^[29]. The high transcriptional level of MYB44 in CK was in an agreement with the lowest level of ABA in CK. MYC2 regulates the termination of JA signaling in tomato^[30]. The transcript of GH3.1 encodes an IAA-amido synthetase that conjugates IAA to amino acids^[31]. In the pathway of plant hormone signal transduction, the expressions of MYC2 (TEA000833.1) and GH3.1 (TEA013731.1) in the pruning-treated samples were much lower than that of CK, which was consistent with the relatively low levels of JA and IAA in CK compared with pruned samples. Circadian rhythms have been associated with shoot architecture in natural settings^[32]. High similarity of HP4.15 and CK was observed in the pathway of circadian rhythm, which are discriminated from other pruning-treated samples. RVE1 (TEA026529.1) was a key discriminator of these two main subgroups. It was reported that RVE1 positively regulated the expression of the auxin biosynthetic gene YUC8 in Arabidopsis, exerting a growth-promoting effect^[33]. The effect of different pruning operations on the circadian rhythm of tea plants can requires further study.

Phytohormones are essential to the growth of tea plants. The growth and transcriptome profiles of tea plants were regulated by different pruning treatments. At the stage of two leaves and one bud, the weight of 100 buds, stem diameter and shoot length of the following spring tea shoots were generally increased due to pruning treatments. Heavy pruning in mid April or early May greatly accelerated the developments of the following spring tea shoots. The levels of phytohormones (IAA, GA1, GA3, and tZ) were significantly increased in the spring buds of tea plants heavily pruned in May. Accordingly, the pathways of plant–pathogen interaction, plant hormone signal transduction and circadian rhythm were enriched under different pruning treatments, based on the results of KEGG. Heavy pruning in April or May, without autumn pruning, caused an early development and elevated weight of 100 buds in the following spring tea shoots. This study provides scientific guidance for plant pruning operations in order to regulate the growth of the following spring tea shoots.

Twenty two-year-old tea plants ('Jinfeng'), grown in the tea plantation of Lanxi Chishan Lake Green Farm Co., Ltd. (Lanxi, Zhejiang Province, China, N 29°11', E 119°48'), were used in this study. Heavy pruning means cutting the main stems back to 40−45 cm from the ground. Deep pruning means cutting the main stems by 10-15 cm from the top. Light pruning means cutting the main stems by 3-5 cm from the top. Different pruning treatments were carried out as follows: (1) Control group (CK), no pruning; (2) Heavy pruning (cutting at the height of 40−45 cm from the ground) on April 15^th 2020, which was termed HP 4.15; (3) Heavy pruning (cutting at the height of 40−45 cm from the ground) on May 8^th 2020, which was termed HP 5.8; (4) Deep pruning (cutting at the height of 10−15 cm from the top) on May 8^th 2020, which was termed DP 5.8; (5) Light pruning (cutting at the height of 3−5 cm from the top) on May 8^th 2020, which was termed LP 5.8; (6) Heavy pruning (cutting at the height of 40−45 cm from the ground) on May 8^th 2020 and light pruning on September 30^th 2020, which was termed HP 5.8+LP 9.30; (7) Heavy pruning (cutting at the height of 40−45 cm from the ground) on May 28^th 2020, which was termed HP 5.28. Field observations were carried out during the spring of the following year. At the bud stage, the buds of different treatments were collected on March 3^rd 2021 for transcriptome, phytohormone and catechin analyses. Three independent biological replicates were collected, and 3–5 tea plants were used for each biological repeat. The buds were immediately soaked in liquid nitrogen for 30 min. The buds belonging to the same biological replicate were mixed on ice and divided into three portions for different analyses. All samples were stored at –80 °C prior to transcriptomic and chemical analyses.

Growth status measurements

Stem diameter was measured in the middle of the third internode from the tip of young shoots, using a vernier caliper.

Shoot length was measured from the tip of the young shoot to the bottom of the young shoot, using a vernier caliper.

Weight of 100 buds (two leaves and one bud): young shoots were plucked and every 100 young shoots were weighed.

Determination of catechin compounds

Tea extract was prepared and analyzed according to our previous method^[34]. Briefly, 0.15 g of the ground tea sample was extracted using 25 mL of 50% (v/v) ethanol solution at 100 rpm and 70 °С for 30 min. After centrifugation (12,000 rpm, 4 °С, 15 min), the supernatants were analyzed by ultra-high-performance liquid chromatography–diode array detector–tandem mass spectrometry (UHPLC–DAD–MS). The conditions of UHPLC (Waters Corporation, Milford, MA, USA) were according to our previously published work^[5]. Catechins were quantified by the authentic standards purchased from Sigma-Aldrich (Shanghai, China).

Transcriptomic and bioinformatic analyses

RNA isolation and sequencing were carried out by Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China) according to the reported method by Ye et al.^[5]. In brief, the total RNA was extracted, enriched and fragmented. The obtained short fragments were reverse-transcripted into cDNA, and second-strand cDNA were synthesized. After purification, end repair, and A addition, the cDNA was ligated to Illumina sequencing adapters. The ligated products were selected by agarose gel electrophoresis, PCR amplified, and purified with AMPure XP beads to obtain the library, and then sequencing was performed on Illumina HiSeqTM 2500. The high-quality clean reads were obtained from filtering raw reads by fastp (version 0.18.0)^[35] through removing adaptor, duplication, and ambiguous sequences (reads with above 10% 'N' rate), as well as low quality reads containing more than 50% of low quality (Q-value ≤ 20) bases, and further were mapped to a reference tea genome of 'Shuchazao'^[36] using HISAT2.2.4, with '-rna-strandness RF' and other parameters set as a default^[36,37]. Mismatches were allowed using the default parameters. Transcript reconstruction was performed on software StringTie v1.3.1^[37]. Fragments Per Kilobase of transcript per Million mapped reads (FPKM) method was used to normalize a gene expression level, using StringTie v1.3.1 software^[37]. The DEGs in the RNA-seq dataset were identified by DESeq2 (ǀlog₂ fold changeǀ > 1, FDR < 0.05) based on read counts^[38]. Significantly enriched KEGG pathways in all the genes compared to the genome background were defined by a hypergeometric test.

RT-qPCR

First-strand cDNA was synthesized from RNA samples (1 µg), using PrimeScript™ RT reagent Kit with gDNA Eraser (TaKaRa Biotechnology Co., Ltd., Dalian, China). Supplemental Table S2 shows the specific primers of selected genes, which was designed by NCBI Primer-BLAST according to the genome sequences of Camellia sinensis 'Shuchazao'^[36]. The qPCR cycling was carried out by Applied Biosystems™ StepOnePlus™ Real-Time PCR System (Applied Biosystems™ ABI, Carlsbad, CA, USA) based on the introduction of PowerUp™ SYBR™ Green Master Mix (Thermo Fisher Scientific Inc., Carlsbad, CA, USA) according to Ye et al.^[5].

Phytohormone analysis

Phytohormone extraction and analysis were carried out by a commercial service company (Gene de novo Biotechnology, Guangzhou, China). The extraction and analysis methods for phytohormones were described in our previous study^[5]. All of the phytohormones were quantified by calculating the area of each individual peak, using authentic phytohormone standards. The HPLC grade of phytohormones, including ABA (≥ 98%), GA1 (95%), GA3 (98%), IAA (≥ 98%), tZ (≥ 95%), JA (95%), tZR (≥ 98%) and melatonin (≥ 95%) were all purchased from Shanghai Yuanye Bio-Technology Co., Ltd. (Shanghai, China).

Data statistics

All the tests were repeated three times and mean value ± SD was presented. The significant difference analysis was performed on the SAS System for Windows version 8.1 (SAS Institute Inc., Cary, NC, USA), using a Tukey test. The Origin 9.1 (Originlab Corporation, Northampton, MA, USA) was used for plotting box plots. Principal component analysis (PCA) based on a correlation matrix was conducted by Minitab 17 statistical software (Minitab. LLC, State College, PA, USA). The heatmap was plotted using z-score values of transcriptomic dataset and drafted on an online platform of OmicShare tools (www.omicshare.com/tools).

Data availability statement

The RNA-seq raw data (Accession number: CRA006814) were up-loaded to BIG data center (https://ngdc.cncb.ac.cn, accessed on 25^th April 2022) under project No.PRJCA009415.