Use of Comparative Transcriptomic Analysis of Different Potato Varieties to Elucidate the Molecular Mechanism Underlying Differences in Cold Resistance


 Background

Among the various abiotic stresses, cold is an essential factor that limits crop productivity worldwide. Low temperature affects the growth, development and distribution of agronomic species around the world. To improve the understanding of the physiological and genetic properties and functions affecting potato cold tolerance, in this study, transcriptomic analysis was performed on two potato strains (HZ88 and LS6) with different cold tolerances that were treated at low temperature for 0, 1, 3, and 6 hours.
Results

Transcriptomic analysis showed that there were large differences between HZ88 and LS6 regarding the expression levels of low-temperature response genes. Notably, HZ88 responds to low-temperature stress, its low-temperature response genes are primarily enriched in plant hormone signal transduction; cutin, suberine and wax biosynthesis; and photosynthesis-antenna proteins. Conversely, the most significant low-temperature response genes of the LS6 strain were determined to be enriched in plant-pathogen interactions, zeatin biosynthesis, and plant hormone signal transduction. The cuticle, composed of a horny waxy layer, is an important protective barrier formed by plants to resist biotic/abiotic stress during the long-term ecological adaptation process, and the HZ88 strain may strengthen its cold resistance by enhancing this physical defence measure. In the LS6 strain, potatoes tend to cope with cold stress by strengthening their immune system and regulating hormone signal transduction. In addition, hormone pathway-related genes (such as ABA), ICE-CBF signalling pathway-related genes, and genes encoding TFs all exhibited different expression patterns between HZ88 and LS6.
Conclusions

To the best of our knowledge, this study is the first to elucidate the genetic mechanisms underlying the difference in cold resistance between the strongly cold-tolerant variety LS6 and the weakly cold-tolerant variety HZ88, thereby establishing a foundation for further analysis and genetic breeding.


Conclusions
To the best of our knowledge, this study is the rst to elucidate the genetic mechanisms underlying the difference in cold resistance between the strongly cold-tolerant variety LS6 and the weakly cold-tolerant variety HZ88, thereby establishing a foundation for further analysis and genetic breeding.

Background
Potato is the world's fourth most important food crop after wheat, maize, and rice [1]. Similar to many other crops [2][3][4][5][6], temperature places strong limitations on the growth and development of potatoes. Low temperature is one of the most common environmental stress factors; it affects potato planting and seedling growth, slows the growth rate, delays tuber maturity, and ultimately reduces rice yield. Cold stress is the main factor leading to the reduction of potato yield in temperate and high-altitude regions [7,8]. Worldwide, approximately 15 million hectares of land are prone to low temperatures, especially in Japan, South Korea, Northeast China and Southwest China [9]. The suitable temperature for potatoes is between 10-21 ℃. At this temperature, potatoes grow rapidly. Beyond this temperature, the growth begins to slow, and the stems and leaves stop growing at 42 degrees. The most suitable temperature for potato owering is 15-17 ℃, and pollen fertility is optimal at this temperature. If the temperature is below 5 ℃ or above 38 ℃, pollen fertility is abolished [10,11]. Therefore, cold stress at the planting and reproductive stages has a negative impact on potato quality and yield. Enhancing the cold tolerance of potatoes is an important goal of potato breeding. This step is necessary not only to maintain potato production in cool areas but also to expand potato planting to northern areas or low-temperature and high-altitude areas. Therefore, the development of hardy potato varieties may bene t food production and contribute to food security and sustainable development.
Plants have evolved complex mechanisms to withstand cold stress, one of which is cold adaptation; in other words, plants have gradually gained freezing tolerance through exposure to nonlethal low temperatures [12]. In this process, a series of comprehensive physiological and biochemical events occur.
At the physiological level, many substances or protective proteins are synthesized in plants, such as soluble sugars, proteases and cold-resistant proteins. These substances are involved in regulating osmotic potential, ice crystal formation, cell membrane stability and reactive oxygen species (ROS) clearance. In the past two decades, many components of the cold stress signalling pathway have been identi ed, including messenger molecules, protein kinases and phosphates, and transcription factors (TFs). The most notable component of this response is the CBF/COR signal pathway [13,14]. The Crepeat binding factor/dehydration response element binding protein 1 (CBF/DREB1) gene is rapidly induced under cold stress and plays a key role in the cold adaptation of plants [15,16]. COR refers to a class of genes regulated by cold stress, such as cold regulation (COR), low temperature induction (LTI) and cold noninducible (KIN), a number of which encode osmotic proteins and cryoprotective proteins to protect plants from frostbite [17][18][19].
Cold tolerance is a very complex quantitative feature that is controlled by many genes. The literature shows that the molecular mechanisms governing cold stress responses have been extensively studied in Arabidopsis. Cold-responsive C-repeat-binding factor (CBF) TFs, including CBF1, CBF2, and CBF3 (also known as DREB1b, DREB1c, and DREB1a), play central roles in plant cold acclimation [20,21]. CBF proteins recognize and bind to the C-repeat/dehydration responsive element (CRT/DRE) motif in the promoters of many cold-responsive (COR) genes, such as RD29a and COR15a, and promote the transcription of these downstream genes [22,23]. The expression of the CBF gene is rapidly induced by cold stress and is both positively and negatively controlled by various TFs. Inducer of CBF expression (ICE1), a bHLH transcription factor, is the best characteristic transcription activator of the CBF gene identi ed to date [24,25]. ICE1 activates CBF by directly binding the CBF gene to the promoter of the promoter under cold stress.
In addition, as a chemical signal, the plant hormone ABA plays a vital role in addressing biological and nonbiological stresses (such as cold) [26]. ABA stimulates several changes in plant physiological, molecular and developmental progressions, resulting in plant adaptation to the stress environment, and it is characterized as an intracellular messenger that has been shown to perform the critical function at the heart of the signalling operation [27]. Abiotic stresses induce the synthesis of ABA which, in turn, activates the expression of stress-related genes and stomatal closure [28][29][30]. In addition, salicylic acid (SA), gibberellin (GA), plant phytohormone JA, and other hormones, such as ethylene and cytokinin (CK), play important roles in cold tolerance. In brief, the levels of plant hormones change under abiotic stress, thereby reducing plant growth.
Although the cold signalling pathway in plants has been extensively studied, the molecular mechanisms underlying cold signal sensing and transduction in plants have not been thoroughly elucidated to date. Especially in potatoes, research on cold resistance mechanisms remains very scarce. In this study, we utilized comparative transcriptomic analysis to characterize the differing responses of the high coldtolerance variety LS6 and the low-hardness variety HZ88 under low-temperature stress. Transcriptomic analysis using RNA-seq was performed for these strains under four kinds of cold treatments at different times points to analyse the difference in cold response between HZ88 and LS6. Many DEGs were identi ed over different times of cold stress, and plant hormone signal transduction was the most recognized molecular pathway. The systematic description of the dynamic changes in the transcriptomic analysis of HZ88 and LS6 under cold stress provided an opportunity to elucidate the molecular mechanisms governing the differences in potato cold resistance.

Results
Characteristics and statistical analysis of de novo transcriptome data Through high-throughput sequencing, a total of 85.05 Gb of high-quality data was generated (Table S1). The raw data of sequencing can be obtained at NCBI through the Bioproject number of PRJNA702083. We rst utilized the published potato genome as a reference for mapping, but the mapping rate was less than 70%. Therefore, considering the incompleteness of the potato genome, we performed transcriptome assembly on the transcriptome data by de novo assembly. A total of 255,612 transcripts were identi ed, of which 138,557 unigenes were further identi ed and predicted by CDS prediction (Fig.S1A). The results showed that 114,422 unigenes were annotated in the NT database, 60,789 unigenes were annotated in the NR database, and a total of 121,892 (87.97%) unigenes could be annotated by known databases (Fig.S1B). In addition, 16665 unigenes could not be annotated in any database, which may represent a further improvement of the potato genome annotation (Fig.S1).

Identi cation of differentially expressed genes
Differentially expressed genes (DEGs) were examined to analyse the genes that may be involved in differences in cold tolerance between distinct potato strains. A total of 1747 DEGs in HZ88 strains were identi ed through HCK vs H1 comparisons, while only 60 DEGs were identi ed between LCK and L1. Among the 1747 DEGs identi ed in the HCK vs H1 comparisons, 1239 downregulated genes and 750 upregulated genes were included ( Fig.1A-B). Among the 74 DEGs in LCK and L1, there were 51 upregulated genes and 23 downregulated genes. Similarly, we employed the same method to identify the different genes at different time points in the HZ88 and LS6 lines. A total of 1001, 1217, 431, 737 DEGs were found in the alignments of HCK and H3, HCK and H6, LCK and L3, and LCK and L6, respectively. Among 1001 HCK and H3 DEGs, 499 upregulated genes and 502 downregulated genes were identi ed, and among 1,217 HCK and H6 DEGs, 583 upregulated genes and 634 downregulated genes were identi ed ( Fig.1C-D). A total of 259 upregulated genes and 173 downregulated genes were identi ed in the comparison between LCK and L3, and between LCK and L6, there were 394 upregulated genes and 343 downregulated genes. In general, the number of upregulated genes in line LS6 was greater than the number of downregulated genes in response to low-temperature stress ( Fig.1A-B), while HZ88 exhibited the opposite result. Especially at one hour of low-temperature stress, as many as 1239 genes in the HZ88 line were downregulated due to low temperature stress, while only 23 genes in the LS6 line were affected by this stress. Comparing the analysis results between different strains showed that LS6 and HZ88 exhibited large differences in their responses to low-temperature stress (Fig.1E), which was re ected in the number of DEGs between the two strains.

Identifying functionally differentiated pathways through GO and KEGG pathway analysis
To further analyse the functional pathways involved in these DEGs, DEGs were evaluated using GO and KEGG pathway analysis. We performed GO analysis on each group of DEGs in the previous results and compared the differences in cold resistance mechanisms between LS6 and HZ88. In DEGs of HCK vs H1 comparisons, we found that upregulated genes are enriched in the biological process of DNA integration and the molecular functions of serine-type carboxypeptidase activity and carboxypeptidase activity, while downregulated genes are enriched in the biological processes of the oxidation-reduction process and single-organism metabolic process and the molecular function of oxidoreductase activity ( Fig.2A). In contrast to the HCK vs H1 comparisons, the upregulated genes in the LCK vs L1 comparisons were enriched in response to auxin biological processes, while downregulated genes were enriched in the biological processes of regulation of nucleic acid-templated transcription and regulation of RNA biosynthesis (Fig.2B). Similarly, we performed GO enrichment analysis of DEGs between HCK and H3, HCK and H6, LCK and L3, and LCK and L6. The upregulated DEGs between HCK and H3 were mostly enriched in DNA integration, serine-type carboxypeptidase activity, and endoplasmic reticulum lumen, while the downregulated DEGs were mostly enriched in the biological process of oxidation-reduction and the molecular functions of haem binding and oxidoreductase activity (Fig.2C). Unlike HCK and H3, most of the upregulated genes between LCK and L3 were enriched in the biological processes of the response to organic substances, response to hormones, and response to abiotic stimuli, while downregulated genes were enriched in the biological processes of the regulation of nucleic acid-templated transcription and regulation of RNA biosynthesis (Fig.2D). Finally, in the comparison of HCK and H6, we found that upregulated genes were all enriched in the biological process of DNA integration, while downregulated genes were enriched in the biological processes of oxidation-reduction, single-organism metabolic processes, and metabolic processes (Fig.2E). In contrast, in the comparison between LCK and L6, most upregulated genes were enriched in to the biological processes of the response to osmotic stress, response to abiotic stimuli, and response to salt stress, while downregulated genes were enriched in the biological processes of the regulation of nucleic acid-templated transcription and regulation of RNA biosynthesis (Fig.2F).
Expression pro le of plant hormone related genes With decreasing temperature, plant growth vigour decreased, and the endogenous hormones produced by plants changed signi cantly, primarily manifested by a large increase in ABA, while such hormones as IAA and GA decreased signi cantly. On the one hand, ABA induces the synthesis of some new proteins related to stress resistance and enhances stress resistance; on the other hand, ABA can promote the closure of stomata, reduce water loss, maintain water balance in the body, and reduce adversity damage. We identi ed 7 important genes related to potato endogenous hormones. Changes in these genes may affect the changes in the content of potato endogenous hormones and may thus affect cold resistance (Fig.3). We identi ed 6 ZEP genes, 2 NCED genes, 17 AAO genes, 4 YUC genes, 3 TDC genes, 1 TAA1 gene, and 3 IPT genes in potatoes subjected to cold stress (Fig.3). Many of these genes exhibited speci c transcriptome change patterns in the HZ88 and LS6 strains. In the LS6 strain, the expression levels of the 6 ZEP genes did not exhibit signi cant differential expression, while in the HZ88 strain, the expression levels of 2 ZEP genes increased strongly in the H1 stage and showed a signi cant upregulation under low temperature stress. expression. The NCED gene also showed a similar change trend. One NCED gene showed a similar change trend between the two varieties, while the other NCED gene was downregulated in the LS6 strain, and it was expressed in the whole stage of HZ88, exhibiting signi cantly upregulated expression. Among the 17 identi ed AAO genes, in the LS6 strain, 2 genes exhibit upregulated expression, while 1 gene was downregulated. However, no signi cant differential expression was observed in the HZ88 strain. Among the 4 YUC and 3 TDC genes identi ed, we identi ed no signi cant DEGs in the LS6 strain, while 1 YUC and 1 TDC gene were signi cantly upregulated in the HZ88 strain.

Differential expression pro les of TFs among strains
TFs play important regulatory roles in cold tolerance. To further explore the mechanisms underlying differences in cold tolerance between the LS6 and HZ88 strains, we identi ed various TFs that are differentially expressed under low-temperature stress. Fifty-seven and 150 differentially expressed TFs were found in the LS6 and HZ88 strains, respectively. The 57 differentially expressed TFs present in LS6 included 23 WRKY family TFs, 6 NAC family TFs, 7 bHLH family TFs, 2 C2H2 family TFs, 1 TCP family transcription factor, and 19 MYB family TFs, while TFs of other families were not differentially expressed (Fig.4A). The 150 differentially expressed TFs found in HZ88 included 43 WRKY family TFs, 50 MYB family TFs, 23 BHLH family TFs, 20 NAC family TFs, 9 TCP family TFs, 3 C2H2 family TFs, 3 FAR1 family TFs, and 1 bZIP family TF (Fig.4B). From the number of TFs, it can be found that the differentially expressed TFs of HZ88 are considerably more abundant than those of LS6. Most of the TFs in LS6 exhibited downregulated expression, and only 1 TCP, 2 bHLH, two WRKY, and four MYB family TFs showed upregulation. In the HZ88 strain, 54% of TFs showed downregulated expression, while 46% of TFs showed upregulated expression. Forty-seven, 47, and 49 upregulated TFs were found in HCK and H1, HCK and H3, and HCK and H6, respectively, while 63, 61, and 46 downregulated TFs were found in HCK and H1, HCK and H3, and HCK and H6, respectively. H3, HCK and H6 were identi ed in the comparison. From this analysis, it can be concluded that WRKY, MYB, NAC, and bHLH family TFs are the most important TFs involved in potato cold resistance, and they also exhibit completely different changes in expression among different varieties.

Expression pro le of the ICE-CBF-COR signalling pathway in different strains
In the past 20 years of research, it has been largely determined that one of the main pathways for the activation of cold response genes is the ICE1-CBF regulatory pathway. Research on CBFs has indicated that overexpression of CBF genes enhances the resistance of plants to freezing and low temperatures.
The multilevel regulation of ICE1-CBF-COR genes is already considered to be one of the important regulation methods governing plant cold resistance. At present, research on this pathway in potato remains very scarce. To further clarify the important role played by the ICE1-CBF-COR pathway in potato cold resistance, we identi ed 1 DREB1 gene, 2 ICE genes, 3 CBF genes and 6 COR genes (Fig.5). We observed that the expression levels of the three CBF genes increased signi cantly under low-temperature stress. With increasing stress time, the expression of CBFs showed a progressive increase and was highest in the H6 and L6 stages. Crucially, we found that the average expression level of CBF in the LS6 strain was approximately two times that of the HZ88 strain. Among these genes, in the HZ88 strain, the fold changes (FC) of the CBF gene at the H6 stage of the differential gene was greater than 200, while in the LS6 strain, the FC was lower, but the expression level was higher. This result may indicate that the CBF gene is directly related to potato cold resistance, and the expression level increased signi cantly with the further increase of the low-temperature stress time. In addition, the remaining 1 DREB1 gene, 2 ICE genes and 6 COR genes also showed differential expression at different stages under low-temperature stress, but no signi cant difference was found between the two potato lines.
Changes in the expression of potato sacchari cation-related genes among different strains Potato tubers degrade starch under low-temperature stress, which leads to the accumulation of reducing sugars, which has a notable impact on the production and storage of potatoes. To analyse the differences in sacchari cation of the two potato lines at low temperature, we identi ed 7 genes that are closely related to low-temperature sacchari cation of potato to explore the differences in low-temperature sacchari cation between the two lines. We identi ed 107 genes related to low-temperature sacchari cation, including 18 sucrose synthase genes, 30 invertase genes, 13 starch synthase genes, 40 amylase genes, 1 APGase gene, and 1 SPS gene, as well as 4 INH genes (Fig.6). We found that the expression level of amylase in the HZ88 line was signi cantly higher than that in the LS6 line, while the expression level of starch synthase in the LS6 line was slightly higher than that in the HZ88 line. The expression level of invertase was higher in the HZ88 strain than in the LS6 strain, while there was no signi cant difference in the expression levels of sucrose synthase. Based on the expression of key genes, it appears that the LS6 strain may have better starch storage capacity in response to low-temperature stress, while HZ88 needs to decompose starch more extensively to cope with cold stress.

Discussion
China is the country with the largest total potato production worldwide [31]. Potatoes are widely consumed because they can provide the human body with abundant calories and are rich in protein, amino acids, vitamins and minerals. Potatoes need appropriate temperatures at all stages of growth, and low temperature is one of the most important limiting factors affecting potato production. In recent years, research on plant cold resistance has made remarkable progress, but the molecular mechanisms governing the responses of potato to low-temperature stress have not been elucidated to date [32]. In this study, we used transcriptomic analysis to determine the differences between the high cold-hardiness line LS6 and the low cold-hardiness line HZ88 in response to low-temperature stress. Speci cally, our analysis showed that the LS6 strain exhibits a weaker expression pattern in hormone-related genes, TFs, ICE-CBFrelated genes and low-temperature glycation-related genes, thereby obtaining better cold resistance than HZ88.
Plant hormone levels are altered under abiotic stresses, and the cause plant growth to decrease. The cold response is one of the most important environmental stresses that induces growth repression and reduces yield. Plant hormone activities induce cascades of cold-responsive genes [32,33]. The plant hormone ABA plays an essential role as a chemical signal in response to biotic and abiotic stresses, such as salt, cold, drought and wounding [34,35]. ABA stimulates several changes in plant physiological, molecular and developmental progressions, resulting in plant adaptation to the stress environment, and it is characterized as an intracellular messenger that has been suggested to perform critical functions as the core of signalling operations [36,37]. In our study, the expression of the two ZEP genes related to ABA biosynthesis was strongly increased in the H1 stage, while the NCED [38,39] gene was downregulated in the LS6 strain. This result indicates that the genes related to ABA biosynthesis in HZ88 show an upwardregulated expression trend as a whole, which may increase the level of ABA in HZ88. In the HZ88 strain, it was possible to limit the growth and development of potatoes by increasing the overall content of ABA, thereby improving cold resistance. However, the expression of related genes in the LS6 strain did not exhibit obvious changes, which indicated that LS6 was less affected when responding to low-temperature stress and did not need to increase the ABA level to limit the growth and development of potatoes. Similarly, we also found that the YUC and TDC genes were signi cantly upregulated in the HZ88 line. These two genes are related to the biosynthesis of IAA, and they also did not exhibit differential expression in the LS6 strain, which may indicate that IAA also plays important roles in the cold resistance of HZ88.
CBFs have been identi ed in all important eld crops and some vegetable species. CBF genes have been identi ed in a wide range of plants [40][41][42]. Intensive research on the expression of cold-response genes in Arabidopsis resulted in the identi cation of the CBF/DREB1 transcription factor family, which plays an essential role in cold acclimation and freezing stress tolerance in plants [41,43]. CBF/DREB1 is a family of TFs that belongs to the ERF/AP2 transcription family [41,[44][45][46]. The upregulation of COR (coldregulated) genes through the activation of CRT/DRE elements in their promoters by CBF proteins is considered to be a small part of the CBF regulon, since intensive studies, which were performed to determine the expression pro ling of CBF, reported that a total of 109 genes were assigned to the CBF regulon [47,48]. We observed that the expression of 3 CBF genes was signi cantly increased under low-temperature stress. With increasing stress time, the expression of CBFs showed a progressive increase. In addition, the average expression level of CBF in the LS6 strain was approximately two times that of the HZ88 strain. Highly expressed CBF improves the cold tolerance of potatoes by activating the expression of downstream related genes. The accumulation of CBF genes in potatoes exhibits a signi cant positive correlation with the time of cold stress, but no such nding has been made in other species.
In addition to the CBF transcription factor family, such TFs as MYB, NAC, bZIP, and WRKY also play important roles in plant cold resistance [49][50][51][52][53][54][55][56][57][58][59][60][61]. The expression of WRKY TFs is induced by external environmental stimuli. Among these genes, 17 WRKY family transcription factor members are induced by low temperature, and most of them are early-response genes. Many WRKY TFs that may play important roles in cold resistance have been identi ed, but their mechanism by which they affect the process of cold resistance has not been elucidated. Researchers cloned 13 WRKY-related genes from the cDNA library of rice leaves, among which the expression of 10 WRKY-related genes were affected by low temperature and other abiotic stress factors. Wang et al. isolated the WRKY46 gene from Chinese cabbage and constitutively expressed the bcWRKY gene through the CaMV 35S promoter in tobacco to improve tobacco resistance to low-temperature stress. In our study, we identi ed 43 and 23 WRKY TFs in the HZ88 and LS6 strains, respectively. In the LS6 strain, these WRKY family TFs largely exhibited signi cant upregulation under cold stress, while in HZ88, these genes were mostly downregulated. This result shows that WRKY, as an important transcription factor for potato cold tolerance, plays an important role in the LS6 strain. NAC family TFs are considered to play important roles in biotic and abiotic stress responses, and some of these genes are considered to be target genes in plant stress engineering. The expression levels of many Brassica napus BrNAC genes and sugarcane SsNAC23 genes are increased under cold, dehydration and pest stresses. The rice OsNAC2 gene can be induced by low-temperature stress, and its overexpression can signi cantly improve the plant's resistance to low temperature and is sensitive to ABA. Similarly, among the NAC family TFs that we identi ed, all genes were upregulated in the LS6 strain, while most genes were downregulated in the HZ88 strain. This result indicates that NAC may play an important role in enhancing the cold resistance of LS6.

Conclusion
Low temperature is one of the serious causal factors of crop yield decline, and potatoes are often subjected to cold stress due to their extensive planting area. In this study, two potato varieties with different cold resistance abilities, LS6 and HZ88, were used to construct a transcriptome map of 1, 3, and 6 hours of low temperature stress. The analysis revealed the differences in low temperature tolerance strategies between the two varieties of potato, and helped us understand the mechanism of yield loss caused by low temperature. The conclusion is that the speci c expression patterns of LS6 and HZ88 in plant hormones, ICE-CBF signaling pathway, WRKY and NAC family TFs, and low-temperature glycationrelated genes lead to differences in their cold tolerance. Speci cally, HZ88 reduces ABA content by reducing the expression of ABA-related genes to limit plant growth to cope with low temperature stress, while LS6 is not affected. WRKY and NAC family TFs were signi cantly up-regulated in LS6 and downregulated in HZ88, which made the WRKY and NAC family TFs we identi ed as possible enhancement factors for low temperature tolerance. At the same time, the expression of genes related to the ICE-CBF signaling pathway is much higher in LS6 than in HZ88, but the ICE-CBF signaling pathway is a key low temperature regulator in both varieties. However, the amylase gene showed a higher expression level in HZ88, which may indicate that HZ88 needs to decompose starch to cope with low temperature stress.
These ndings provide a variety-speci c gene expression network for potato under low temperature stress, which helps to understand the mechanism of low temperature stress response, so as to reduce the reduction of potato yield caused by low temperature, and further guide potato breeding.

Plant Materials and RNA quanti cation
In this study, the strongly cold-tolerant strain (LS6) and the weakly cold-tolerant strain (HZ88) were employed to elucidate the molecular mechanism governing the differences in low-temperature tolerance between potato strains in the short term. Each potato line was grown under 4 conditions: low temperature at 0 (CK) and 1, 3 and 6 hours. Each time node employed 3 biological replicates, and nally, a total of 24 samples were used for RNA-seq analysis. RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was checked using a NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using a Qubit® RNA Assay Kit in a Qubit® 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Library preparation and transcriptome sequencing
Each sample contained 1.5 µg RNA as the preparation material for RNA-seq, and construction of the sequencing library strictly followed the recommendations of the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA). Through puri cation, reverse transcription and ampli cation of total RNA, a primary cDNA library was constructed. Next, the AMPure XP system (Beckman Coulter, Beverly, USA) was used to screen cDNA fragments with a length of 150~200 bp, and PCR was performed with Phusion High-Fidelity DNA polymerase, universal PCR primers and Index (X) Primer. Finally, PCR products were puri ed (AMPure XP system), and library quality was assessed on the Agilent Bioanalyzer 2100 system to prepare for sequencing.
Clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq platform, and paired-end reads were generated.

Transcriptome assembly and expression quanti cation
The left les (read1 les) from all libraries/samples were pooled into one large left.fq le, and the right les (read2 les) were pooled into one large right.fq le. Transcriptome assembly was accomplished based on the left.fq and right.fq using Trinity with min_kmer_cov set to 2 by default, and all other parameters were set to default.

Differential expression analysis and enrichment analysis
Differential expression analysis of two conditions/groups was performed using the DESeq R package (1.10.1). DESeq provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P-values were adjusted using Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with adjusted P-values were subjected to Gene Ontology (GO) enrichment analysis of the DEGs, which was implemented by the GOseq R package-based Wallenius noncentral hypergeometric distribution, which can adjust for gene length bias in DEGs. KEGG (Kyoto Encyclopedia of Genes and Genomes) is a database resource for characterizing the high-level functions and utilities of biological systems, such as cells, organisms and ecosystems, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/). We used KOBAS software to test the statistical enrichment of differentially expressed genes in KEGG pathways. Availability of data and materials All data and materials are presented in the main paper and additional supporting le.

Ethics approval and consent to participate
There are no ethical issues involved.

Consent for publication
Not applicable.