Transcriptional reprogramming of the bud-mutation loquat YongLu emerging freezing resistant

Background: Freezing seriously affects loquat, an originally subtropical fruit. Here, a wide-spread cultivar loquat bud mutation ‘Yonglu’ (YL) was found to be more resistant to freezing than the parental ‘Ruantiaobaisha’ (RT), especially the fruits. The freezing resistance mechanism of YL was analyzed by physiological and transcriptomic methods. Results: After freezing treatment, only the fruits of YL showed a signicant accumulation of proline (which increased by 1.5-fold). Interestingly, abscisic acid (ABA) accumulation in YL showed no difference after freezing treatment. Moreover, the stomatal density, area, and apertures of YL signicantly decreased. All these results suggested that the proline content and stomatal closure contributed to the enhanced freezing tolerance. Further, the transcriptional proles of freezing tissues were obtained and compared. Functional enrichment analysis showed that aside from basic nutrients and energy metabolisms, such as citrate cycle and oxidative phosphorylation, the fruits and owers of YL have different systems to enhance freezing tolerance, suggesting a progressive freezing-resistant network. Signaling transduction of nearly all hormones, including ABA, was signicantly differentially expressed. Several key genes, such as PP2C and SnRK2, which are crucial for transducing ABA signals, were signicantly down- and upregulated, respectively. But the genes of neither ABA biosynthesis nor ABA transduction pathway signicantly differentially expressed between two freezing varieties, suggesting they contribute limit in enhancing freezing tolerance of YL. Conclusion: Concludingly, the bud mutation YL loquat has stronger freezing tolerance than its parental RT due to the increasing proline and more effective regulation of stomatal closure. Furthermore, the bud mutation YL loquat reprogrammed its entire transcriptional proles to enhance their abilities to adapt to freezing stress. However, ABA biosynthesis, ABA transduction pathway and CBF pathway, which are common reasons in cold response, were not the main reason for the increasing freezing tolerance of YL. There must be unknown regulation mechanisms for the bud mutation to freezing resistant.


Introduction
Loquat (Eriobotrya japonica Lindl.) is an evergreen fruit tree originating in southern China, and China accounts for more than 80% of the world's planting area and fruit production of loquat [1][2][3]. Loquat blooms in winter. Its ower buds can tolerate a low temperature of − 6 °C. However, the ovule of its fruit becomes yellow and sheds at − 3 °C, resulting in heavy losses [4]. Therefore, breeding freezing-resistant varieties is important for the industrial development of loquat [1,5,6].
To reduce the damage of low temperature, plants form a series of molecular and physiological response pathways, such as the accumulation of compatible osmolytes (soluble sugars, glycine betaine, and proline), antioxidants, soluble proteins, and a series of freezing response regulatory proteins [7][8][9][10]. For example, C-Repeat/Drought-Response-Element binding factor (CBF/DREB) transcription factors were considered as an important regulator in responding to freezing and cold stress [11][12][13][14]. In Poncirus and Citrus, the expression (time and number) of CBF/DREB is closely related to the tree varieties with large differences in freezing resistance. The correlation indicates that the CBF pathway plays an important role in freezing resistance [15]. In previous studies on the freezing resistance of fruit trees, some cold-resistant proteins were found, and their functions were clari ed [16,17]. For example, dehydrins (DHNs) are a family of plant proteins that are typically induced in response to stress conditions, such as low temperatures, high salinity, and drought, causing cellular dehydration [18]. By comparing the transcriptome data of freezing-sensitive 'Ninghaibai' (FS-NHB) and freezingtolerant 'Jiajiao' (FT-JJ), Eriobotrya japonica DHN proteins (EjDHNs) were found to be closely related to freezing tolerance; these proteins might participate in the cryoprotection of the plasma membrane, preventing its dehydration in response to freezing damage [19]. In another study, EjDHN1 was overexpressed in tobacco, and the transgenic plants showed increased tolerance to cold stress; EjDHN1 could protect cells against oxidative damage caused by ROS production under cold stress [20].
All these genes are regulated by several plant hormones. Plant hormones are low-molecular-weight molecules that are the major responding signals when plants sense abiotic and biotic stresses, including cold and freezing; they are also involved in growth, development, and stress response [21][22][23]. Cold/freezing stress changes the content of most hormones in plants [23]. Interestingly, cold/freezing stress represses the contents of hormones involved in growth, such as gibberellin (GA), brassinosteroid (BR), and cytokinin (CK), and induces those of defense hormones, including abscisic acid (ABA), jasmonic acid (JA), and salicylic acid (SA) [23]. Thereafter, the hormone signaling transduction pathways are triggered along with the changes of hormone contents. At present, several hormone signaling transduction pathways have been identi ed. For example, BR, CK, and ABA play positive roles in enhancing cold/freezing tolerance of plants [23]. Interestingly, although the accumulation of BR and CK is repressed during cold stress, they can promote cold tolerance [24,25]. Among hormones, ABA plays central roles in enhancing freezing tolerance [26,27]. Moreover, GA plays opposite roles to ABA in cold response as a low content of GA would promote freezing tolerance in plants [28].
However, studies on the cold resistance mechanism of loquat are limited due to the scarcity of cold-resistant varieties.
Thus, it is di cult to clarify the molecular regulation mechanism of freezing resistance, which hinders the selection of key genes for freezing resistance [29,30]. In the present study, we selected the cold-resistant variety 'Yonglu' (YL) through bud transformation from the wild-type loquat variety 'Ruantiaobaisha' (RT). YL loquat was found and quali ed in 2007, and it has been widespread cultivated for over 5,000 acres in Zhejiang Province since 2018. After years observations, YL has series excellent agronomic traits, such as high fruit eld and great taste. Furthermore, the freezing resistance of its fruits was signi cantly stronger than that of RT and other varieties. We conducted RNA sequencing to detect physiological indicators and selected the freezing-resistant genes.

Results
YL fruits and owers had less freezing rate under low temperature YL is a bud mutation found from a 35-year-old RT loquat in Tongyu Street, Bututang Village, Luqiao District, Taizhou City, Zhejiang Province in 2003. After years of observation, signi cant differences in fruit shape, seed number, and edible rate were found between YL and its mother plant RT. For example, the number of seeds per fruit in YL is 1.03 less than that in RT ( Figure S1).
Given that loquat is originally a subtropical fruit tree, freezing tolerance became the major factor of fruit production. Through years of investigations on the annual occurrence of freezing damage, bud mutation YL exhibited strong freezing tolerance than its parent RT, especially the fruits ( Fig. 1a and 1c).
Under low-temperature conditions, RT owers were susceptible to freezing damage, which was characterized by deformed owers, incomplete differentiation of ower buds, and local browning. Thus, according to these abnormal phenotypes, the freezing owers of RT and YL were counted. Finally, the freezing rates of YL owers were only 1.70% at those years with freezing winter (2010, 2014 and 2017), whereas those of RT were 6.30% (Fig. 1d).
At the fruit growth stage, RT fruits wilted and their color became purple brown after freezing, and these brown fruits were counted as freezing samples. Totally, the freezing rate of YL fruits was 8.60%, which was extremely lower than that of RT fruits (55.20%) after freezing damage, indicating that YL fruits were more resistant to cold than RT fruits ( Fig. 1c and 1e).
In conclusion, YL showed enhanced freezing tolerance than RT. The fruits of YL had the strongest ability against freezing; however, the owers also had a modest freezing tolerance. Proline Was Signi cantly Accumulated In YL Fruits Proline and ABA accumulation is commonly considered as the main response of plants to defend against freezing damage [31]. A study has shown that plants would proactively decrease the content of GA to promote freezing tolerance [28]. To clarify the main factors behind the enhanced freezing tolerance in YL, the contents of proline, ABA, and GA were quanti ed.
For proline, only YL fruits showed a signi cant increase (almost 1.44-fold) after freezing treatment (Fig. 2a). RT fruits showed no difference in proline content after freezing treatment, but the leaves of both YL and RT had an insigni cant reduction in proline content. Thus, proline accumulation might be bene cial for enhancing freezing tolerance in YL fruits.
For ABA and GA, only RT showed signi cant changes. ABA and GA were increased in the leaves of RT, but GA content was decreased in RT fruits (Fig. 2b). In general, increased ABA and decreased GA are bene cial to freezing resistance [27,28]. In YL leaves and fruits, ABA or GA showed no difference.

YL Had Stronger Actions On Stomatal Closure
Plant stomata are very important for the exchange of water and CO 2 ; thus, they are important structures that protect cells from cold [32]. Here, we observed stomatal structures with a scanning electron microscope. After freezing stress, the stomatal density of RT and YL leaves was obviously decreased, and only YL showed a signi cant decrease (Fig. 3).
In fruits, a single stoma and total stomatal area of YL fruits were signi cantly reduced after freezing (Fig. 4). The stomatal area of YL fruits was 1.77 ± 0.12 and 1.16 ± 0.12 µm 2 before and after freezing, respectively. The stomatal apertures of YL fruits were obviously smaller than those of RT fruits as the long and short axes were signi cantly reduced by 0.8-fold. Conversely, RT fruits did not exhibit any signi cant differences.
Scanning electron microscopy results showed that stomatal closure is probably an important adaptation strategy to enhance freezing tolerance in YL fruits.
Transcriptome analysis revealed that the fruits and owers of YL applied different freezing resistance systems In order to clarify the freezing tolerance mechanism of YL, especially its fruits, three freezing tissues were subjected to high-throughput sequencing in order to screen the differentially expressed genes (DEGs) of YL and RT. A total of 268,698 unigenes were obtained with an N50 of 1100 bp, and they were annotated against NR, NT, Uniprot, Kyoto Encyclopedia of Genes and Genomes (KEGG), eggNOG, and Gene Ontology (GO) databases. Their expression levels were calculated by the reads per kilobase million mapped reads (RPKM) method ( Figure S2, Tables S1 and S2). Interestingly, the transcription pro les of owers were quite different between YL and RT, whereas the fruits and leaves exhibited similar expression patterns, suggesting that these differences were caused by the freezing tolerance system ( Figure S2e).
Given that all three tissues of YL exhibited strong freezing tolerance to different degrees (fruits > owers > leaves), integrative analysis of DEGs was helpful to clarify the antifreeze mechanism of YL, especially the fruits. DEGs were selected by DEGSeq2 software (R package, https://bioconductor.org/packages/release/bioc/html/DESeq2.html) with a threshold of |log 2 ratio(fold change)|≥1 and q < 0.05 (FDR corrected) (Fig. 5, Tables S3-S5). Finally, a total of 807 DEGs were obtained in all tissues from YL and RT, suggesting that they play key roles in freezing tolerance (Fig. 5a). To verify the transcriptome data, 20 genes were selected to test their expression levels in all samples by qRT-PCR. The results were highly consistent with the transcriptome data (Fig. 5d).
According to the enrichment analysis of functional annotation, freezing leaves only enriched "Ribosome" (Fig. 5b). However, aside from "Ribosome", owers also mostly applied series basic nutrients and energy metabolisms ("Citrate cycle," "Glycolysis," amino acid metabolisms, and "Oxidative phosphorylation") and "Phagosome," suggesting that the owers were able to respond to freezing stress by alternating some pathways of growth and development when the stress was not severe for them ( Fig. 5b and 5c).
The situation in the fruits was quite different (Figs. 5 and 6, Table S3). Aside from several crucial nutrient and energy metabolism pathways including the "Citrate cycle," "Glycolysis," and "Oxidative phosphorylation," fruits recruited more signal transduction and defense pathways, such as "Plant-pathogen interaction," "Plant hormone signal transduction," and "MAPK signaling pathways-plant" (Figs. 5b and 6d). Most enrichment pathways were downregulated, but approximately 18.13% (33) of signi cantly different GO terms (182) were upregulated in YL fruits, such as "DNA integration" and "RNA-direct DNA polymerase activity" (Fig. 6c, Tables S6-S8). Moreover, the rst 10 upregulated GO terms were mainly related to DNA synthesis and modi cation.
Taken together, these results showed that different tissues of YL applied different freezing resistant systems due to unknown reasons. This might explain why different phenotypes were observed in different tissues after freezing treatment.
Neither ABA biosynthesis nor ABA pathway might be the main reason for the freezing resistant According to our results, the main antifreeze hormone, ABA, was not signi cantly accumulated (Fig. 2). Although the "Plant hormone transduction" pathway including the ABA signaling pathway signi cantly differentially expressed and enriched, ABA biosynthesis and signaling transduction pathway seemed not to be the main reason for freezing resistant of YL fruits (Figs. 5 and 6).
According to RNA-seq data, NCED (9-cis-epoxycarotenoid dioxygenase), a key biosynthesis gene of ABA, expressed no difference in freezing YL fruits than RT fruits (Table 1). Since the content of ABA showed no difference, suggesting that accumulation of ABA was not prior to enhance freezing resistant. Furthermore, genes of ABA pathway were investigated. Only two genes, PP2C (protein phosphatase 2C) and SnRK2 (SNF1-related protein kinase 2), which play decisive roles in ABA signal transduction, were all differentially expressed (Table 1). PP2C is a negative regulator in ABA signal transduction [33]. A total of 15 PP2C genes were signi cantly differentially expressed, and they were all downregulated in freezing YL fruits compared with freezing RT fruits (Table 1). Contrary to PP2C, SnRK2 is a positive regulator but is inhibited by PP2C in ABA signal transduction [33]. According to the transcriptome data, only one SnRK2-like (TRINITY_DN83302_c0_g2, closest to SnRK2.8) was signi cantly differentially expressed, and it was upregulated by 2-fold in YL fruits exposed to cold stress (Table 1). But SnRK2.6 was the main regulator to promote stomatal closure by increasing the phosphorylation of SLAC1 (slow anion channel 1) and repressing the phosphorylation of KAT1 (potassium channel protein) [34]. In YL fruits, the expression levels of representative SLAC1 (TRINITY_DN80888_c1_g1, TRINITY_DN84572_c1_g4) and KAT1 genes were changed to some degree but not signi cant (Table 1).
Moreover, we investigated the marker genes of ABA pathway and CBF pathway, and only CBF4 (TRINITY_DN84340_c0_g1) down-regulated for 0.44-times signi cantly in YL fruits, whereas there was no difference of Rd29b (TRINITY_DN83533_c2_g1 and TRINITY_DN73189_c0_g1), KIN1 (TRINITY_DN85365_c0_g1) and LEA76 (TRINITY_DN75564_c0_g1 and TRINITY_DN61955_c0_g1) between two freezing fruits (Table 1). Due to genomic loss of loquat, RAB19 and COR15B orthologs were not found in these samples. Thus, it seemed that ABA pathway was not the regulation system of YL fruits to promote stomatal closure.
Above all, neither ABA biosynthesis nor ABA pathway might be the main reason for the freezing resistant.

Discussion
Freezing stress is an inhibitory factor for the growth and development of fruit trees, especially subtropical plants such as loquat. Freezing damage would cause heavy loss of fruit production, leading to economic losses of local farmers. Thus, breeding loquat varieties with strong freezing tolerance has important economic values [1,5].
Bud mutation is an effective and commonly used technology to breed new varieties of fruit trees. In 2003, we found a bud mutation named YL from RT, which is one of the best quality cultivated loquats in Southern China. After years of observation, we found that YL has greater economic values than RT due to its excellent traits, such as heavier fruits, less seeds, and lower cracking rate. In particular, YL has a stronger freezing tolerance than parental RT (Fig. 1). In 2010, 2014, and 2017, severe frosts occurred in local nurseries, and the damage of the two varieties was evaluated. According to the results, the freezing rate of YL fruits was 8.6% whereas that of RT fruits was 55.2% suggesting that bud mutation enhanced the freezing resistance of fruits (Fig. 1e). The freezing rate of owers was reduced from 6.3-1.7%, suggesting that RT owers has a modest freezing tolerance (Fig. 1d). For leaves, no obvious difference was observed according to the phenotype. However, according to the transcriptome data, YL leaves had altered pathways, such as "Ribosome," which might be involved in cold response (Fig. 5b, Tables S7 and S8). Several ribosomal proteins were considered to be involved in responding to cold, and "Ribosome" might be related to freeze resistance [35].
All these results suggested that YL was more tolerant to freezing than RT. Because different tissues showed different degrees of freezing tolerance, the mutations might be bene cial for owers and fruits to recruit antifreeze activities.
To clarify the freezing resistance mechanism of YL, a series of physiological tests was conducted. According to the results, the content of proline was signi cantly altered, and it might be closely related to the antifreeze activities of YL (Fig. 2a).
In most plants, the ABA content is induced by low temperature and then plays central roles in enhancing freezing tolerance in plants [12,23,26,27,34]. However, according to our results, only RT leaves accumulated ABA after cold treatment, but the ABA content in RT fruits was reduced (Fig. 2b). The ABA content increased in YL leaves and fruits, but the difference was not signi cant, suggesting that the ABA content was not the reason why YL has better antifreeze system (Fig. 2b). Also, the key biosynthesis and metabolizing genes, NCED and CYP707A, were not differentially expressed either (Table 1). These results suggested that ABA accumulation is not the main reason why the bud mutation gain freezing resistant.
When plants are exposed to low-temperature stress, the number of plant stomata and stomatal density would decrease [23,36]. Agurla et al. found that the guard cells on the stomata can sense abiotic and biotic pressures, thereby responding quickly to shut down under stress conditions [37]. The results of this study are consistent with those of previous studies. Stomatal closure seemed to be the determining factor that made YL more resistant to freezing (Figs. 3  and 4). RT leaves did not show signi cant damage after cold treatment; however, they exhibited similar stomatal closure activities as YL leaves (Fig. 3). In fruits, only YL fruits exhibited signi cant and obviously reduced stomatal area and smaller stomatal apertures (Fig. 4).
All these results suggested that YL has stronger freezing tolerance because they have more strict regulation systems in stomatal closure.
Stomatal closure contributed to freezing resistance in YL, especially in the fruits. To further clarify the molecular mechanism of YL in response to freezing, transcriptome data were obtained (Tables S1 and S2). According to the results, different tissues exhibited different levels of freezing tolerance when plants respond to low-temperature stress (Fig. 5). Defense systems such as "Plant hormone signal transduction" and "MAPK signaling" were extremely enriched in YL fruits although these pathways were also differentially expressed in YL leaves and owers (Tables S6-S8). "Ribosome" and "Phagesome" were differentially expressed and enriched in the leaves and owers, suggesting that they contributed to freezing resistance. Several studies revealed that ribosomal protein accumulation and autophagy are related to freezing or drought response in plants [35,38]. Thus, the detailed relationship between the two pathways and freezing resistance needs to be elucidated.
Neither the ABA content was not signi cantly increased in YL fruits, nor the key genes involved in ABA signal transduction pathway barely signi cantly differentially expressed after exposing to cold stress (Table 1). Only PP2C and SnRK2 showed slightly differentially expressed, which are considered to be crucial factors in transducing ABA signals and regulating stomatal closure [33,34,37]. However, the differentially expressed SnRK2 in loquat was closest to SnRK2.8 while not SnRK2.6 which was the main factor to regulate stomatal closure. The downstream genes, SLAC1 and KAT1 directly activated and repressed by PP2C-SnRK2 and are responsible for stomatal closure, showed no difference in YL fruits either. Furthermore, the expression patterns of marker genes of ABA pathway (CBF4 and Rd29b) also showed that ABA signal transduction pathway was not enhanced in YL fruits, suggesting that neither ABA accumulation nor ABA pathway is responsible for the freezing tolerance of YL. All these results suggested that there are other unknown molecular mechanisms in stomatal closure and antifreeze activities in YL fruits.
Interestingly, CBF transcription factors were veri ed to be induced by low temperature. The CBF homologues (TRINITY_DN84340_c0_g1 and TRINITY_DN80208_c2_g1) showed a downregulation tendency in YL tissues compared with RT tissues [11,14]. Furthermore, 12 out of 13 EjDHN1-like genes also showed a downregulation tendency in YL fruits (Table 1) [20]. All these changes of gene expression showed that YL has a speci c antifreeze mechanism.
In summary, the bud mutation YL has stronger freezing tolerance than its parental RT. The increased proline and effective regulation of stomatal closure contributed in enhancing freezing tolerance. Interestingly, the accumulation and transduction of ABA molecules, and recruitment of antifreeze regulators, such as CBFs, were not the rst option for YL to enhance freezing tolerance. There must be other unknown pathways contributing to stomatal closure and enhancing freezing tolerance in YL, such as "Ribosome" and "Phagosome." More detailed molecular functional veri cation should be conducted to further elucidate the relationships between these pathways and freezing tolerance.

Materials And Methods
Plant materials and freezing treatment Thirty fruits of YL and RT grown to 30 days were randomly picked and placed in an incubator (Binder kbf720, Germany) to freeze the sampled fruits for 12 h at − 5 °C. Thirty leaves were randomly taken from each plant, and the sampled leaves were subjected to freezing injury treatment (3 °C to 4 °C) for 12 h in the incubator. The above tests were biologically repeated three times.

Investigation Of The Freezing Rate Of Fruits And Flowers
The freezing conditions of RT and YL were investigated. Flower frostbite mainly manifests on the petals. At years of 2010, 2014 and 2017, the winter were freezing, we investigated the freezing rates of owers of RT and YL. RT and YL were planted in same nursery in Taizhou, Zhejiang. The entire brown fruits were selected as the statistical object of the frozen damaged fruits, and the slight browning of the epidermis was not included. Meanwhile, the growth of the fruit stalk was also investigated.
Observation and statistics of the stomatal structures of leaves and fruits A total of 20 stomatal structures of leaves and fruits were observed with a scanning electron microscope (DMI3000, Leica Microsystems, Germany) after cold/freeze treatment. Image-Pro Plus was used to measure the stomatal short axis, long axis, density, and area.

Detection of ABA and GA hormones and proline content in leaves and fruits
The leaf and fruit tissues of RT and YL were accurately weighed before and after freezing, and 9-times PBS was added to make a 10% (weight-volume ratio) tissue homogenate. Then, the homogenate was centrifuged at 2500 rpm for 10 min and sampled for 50 µl. The ABA content of each sample was determined in accordance with the instructions of the Plant Abscisic acid (ABA) ELISA Kit (Cusabio BIOTECH Co. Ltd., Wuhan, China). The GA content of each sample was determined by using the Plant Gibberellic acid (GA) ELISA Kit (Cusabio BIOTECH Co. Ltd., Wuhan, China). The protein concentration was determined by the Coomassie Brilliant Blue method, and the proline content in each sample was determined using the colorimetric method. The above tests were biologically repeated three times.

Total RNA Isolation
The samples were sent to Annoroad Gene Technology (Beijing, China) for high-throughput sequencing. The CTAB method for RNA isolation was conducted in accordance with the method of Shan et al. [39] with slight modi cations. The total RNAs of each leaf, ower, and fruit of YL and RT were extracted. The RNA samples were then tested for completeness by 0.8% agarose gel electrophoresis. The NanoDrop-2000 spectrophotometer (Thermo Scienti c, USA) was used for concentration detection, the Kio K5500 spectrophotometer (Kio, Beijing) was used to test sample purity, and the Agilent 2100 RNA Nano 6000 Assay Kit (Agilent Technologies, CA, USA) was used to detect the integrity and concentration of RNA samples. All treatments and different tissues have three bio-repeats.
CDNA Library Construction, High-throughput Sequencing, And Gene Annotation The mRNA was enriched with magnetic beads with Oligo (dT). Then, the fragment buffer was added to the resulting mRNA to make the fragment shorter. The fragmented mRNA was used as a template to synthesize the rst strand of cDNA with six base random primers. Next, the buffer, dNTPs, RNaseH, and DNA Polymerase I were added to synthesize the second strand of cDNA. After puri cation using the QIAQuick PCR kit (QIAGEN, Germany) and elution with EB buffer, the puri ed double-stranded cDNA was eluted and then subjected to end repair, base A, and sequencing adapter treatment. Then, the target size fragments were recovered by agarose gel electrophoresis. Finally, the PCR ampli cation was conducted to complete the entire library preparation.
The constructed library was sequenced with Illumina HiSeq 2500 platform. The sequencing strategy was pair-end sequencing, and each raw read was 150 bp in length (PE150). The linker and low-quality reads were used in the original data (Q ≤ 19, bases ≤ 15%, or N-containing ≥ 5%; raw reads at either end were low-quality reads; double-end reads were removed) to obtain high-quality clean data. The Trinity software was used to assemble unigenes, and the Trinotate software was used to annotate ORF and unigene. Thereafter, the predicted sequences with six channels including NT, NR, Swiss-Prot, eggNOG [40], KEGG [41,42], and GO were annotated [43].

DEG Classi cation And Annotation
The RPKM method was used to calculate the expression level of each unigene. DESeq2 was used to analyze the differential expression of the two groups of samples with biological duplication [44]. The Wald test was used to calculate the non-differential expression of each gene in the two groups of samples. The q-value (FDR) was obtained by performing multiple test correction on the p-value using the Benjamini-Hochberg (BH) methods. q < 0.05 and |log2FC|>1 were selected as the selection criteria for differential genes [45]. Based on the log2FC value of the gene, K-means cluster analysis was used to visualize the expression pattern of differential genes using MeV (Multiple Experiment Viewer) software [46].
Differential genes were subjected to GO (http://geneontology.org/) and KEGG (http://www.kegg.jp/) enrichment analyses. According to the number of DEGs contained in each GO term, Fisher's exact test was used to nd GO terms or KEGG pathways that were signi cantly enriched to the entire genomic background. After multiple test correction of Pvalues using the BH method, q < 0.05 was selected as the GO Term or KEGG Pathway screening criterion enriched in DEGs.

QRT-PCR
The total RNA of the transcriptome sequencing was used as a template, and reverse transcription was performed using TUREscript 1st Stand cDNA SYNTHESIS Kit (Aidlab, Beijing). The quantitative primers of the selected unigenes were designed using Primer Premier 5.0 software (http://www.premierbiosoft.com/primerdesign/) ( Table S9). The primers are listed in Table S9. Then, 2 × SYBR Green Mix (Beijing unique biotechnology company, Beijing, China) was used in qRT-PCR, and the uorescence was quanti ed with QTOWER2.2 (Analytik Jena, Germany). Gene expression was calculated with β-actin as internal reference, and the 2 −△△Ct method was used to calculate the relative gene expression of each sample. All treatments and different tissues have three bio-repeats.

Declarations Con ict of Interest
All the authors declared that they have no con icts of interest to this work.