Expression of Genes Related to Plant Hormone Signal Transduction in Jerusalem Artichoke ( Helianthus tuberosus L.) Seedlings under Salt Stress

: Background: Jerusalem artichoke ( Helianthus tuberosus L.) is moderately tolerant to salinity stress and has high economic value. The salt tolerance mechanisms of Jerusalem artichoke are still unclear. Especially in the early stage of Jerusalem artichoke exposure to salt stress, gene transcription is likely to undergo large changes. Previous studies have hinted at the importance of temporal expression analysis in plant transcriptome research. Elucidating these changes may be of great signiﬁcance to understanding the salt tolerance mechanisms of it. Results: We obtained high-quality transcriptome from leaves and roots of Jerusalem artichoke exposed to salinity (300 mM NaCl) for 0 h (hour), 6 h, 12 h, 24 h, and 48 h, with 150 and 129 unigenes and 9023 DEGs (differentially expressed genes). The RNA-seq data were clustered into time-dependent groups (nine clusters each in leaves and roots); gene functions were distributed evenly among them. KEGG enrichment analysis showed the genes related to plant hormone signal transduction were enriched in almost all treatment comparisons. Under salt stress, genes belonging to PYL (abscisic acid receptor PYR/PYL family), PP2C (Type 2C protein phosphatases), GH3 (Gretchen Hagen3), ETR (ethylene receptor), EIN2/3 (ethylene-insensitive protein 2/3), JAZ (genes such as jasmonate ZIM-domain gene), and MYC2 (Transcription factor MYC2) had extremely similar expression patterns. The results of qRT-PCR of 12 randomly selected and function known genes conﬁrmed the accuracy of RNA-seq. Conclusions: Under the inﬂuence of high salinity (300 mM) environment, Jerusalem artichoke suffer serious damage in a short period of time. Based on the expression of genes on the time scale, we found that the distribution of gene functions in time is relatively even. Upregulation of the phytohormone signal transduction had a crucial role in the response of Jerusalem artichoke seedlings to salt stress, and the genes of abscisic acid, auxin, ethylene, and jasmonic acid had the most obvious change pattern. Research emphasized the regulatory role of hormones under high salt shocks and provided an explorable direction for the study of plant salt tolerance mechanisms.


Introduction
More than 800 million hectares of land worldwide are affected by salinity, and the utilization of saline soil is a key issue in agriculture [1]. Moreover, saline soils may be inhospitable to plants. Most crops are glycophytes, meaning they are relatively sensitive to salinity [2][3][4]; hence, finding a suitable crop to grow on saline land is difficult.
Jerusalem artichoke (Helianthus tuberosus L.), a plant species with a variety of potential uses, is moderately salt tolerant [5] and is suitable for saline soil cultivation [6][7][8], which has extraordinary significance for carbon neutrality. The current research on Jerusalem artichoke focuses on breeding for greater biomass production [9][10][11], but less work has been done on its salinity tolerance. Different tissues of Jerusalem artichoke have differential salt tolerance, increasing from tubers to stem and leaves. Under salt exposure, no significant difference was found in Cl − content in tubers, stems, and leaves of Jerusalem artichoke, but Na + content decreased from tubers to leaves, suggesting restricted Na + transport into stem and leaves [12].
Different Jerusalem artichoke genotypes may have differential salt tolerance mechanisms. For example, genotype Nanyu No. 1 adapts to NaCl stress by reducing net photosynthetic rate, stomatal conductance, and transpiration while maintaining a high degree of H + -ATPase activity, whereby genotype Qingyu No. 2 accumulates proline and maintains chlorophyll concentration [13]. Although the amount of literature on salt tolerance of Jerusalem artichoke is increasing, the underlying mechanisms remain unclear.
Many studies have explained the mechanisms of plant salinity tolerance based on physiological, biochemical, and molecular characterization [14]. With the development of high-throughput sequencing, the RNA-seq technology has increasingly been used to study transcriptome of many plant species, such as wheat, barley, salvia, soybean, pepper, and other economic crops [15][16][17][18][19][20], but there is relatively little research on Jerusalem artichoke [10,21,22] and even less so regarding salt stress [21]. Therefore, it is important to use RNA-seq to characterize salt tolerance mechanisms of Jerusalem artichoke. Comparative transcriptomic research focuses on time-course experiments to identify the expression patterns of genes. Such time-series research has been used in the exploration of various biological processes [23][24][25][26].
The important roles of plant hormones in resistance to abiotic stresses have been emphasized [27][28][29], including mediation of salt stress to regulate plant adaptation [30]. Many studies have showed that the expression of genes related to phytohormones depends on the severity and/or duration of abiotic stresses [31,32], indicating that changes in the expression of phytohormone-related genes need to be characterized with respect to stress duration. Therefore, the aim of the present study was to use the transcriptome samples of tissues of Jerusalem artichoke exposed to salt stress for different durations to characterize the expression patterns of genes (especially those related to phytohormone signal transduction).

Plant Material Preparation and Salt Treatment
Jerusalem artichoke (Helianthus tuberosus L. cv. Nanyu No. 1 (NY-1)) tubers were collected in Dafeng (Jiangsu Province, China). Each bud eye was sliced into pieces at least 1-cm thick when breeding; the pieces were surface-sterilized and placed in quartz sand (30 cm × 20 cm enamel tray) moistened with distilled water. The initial cultivation was done at 25 • C in the dark using an incubator (NingBo SaiFu, PGX-350B, Ningbo, Jiangsu Province, China). After sprouting, cultivation continued at 25 • C and 12 h light/12 h dark using 1/2 Hoagland solution to moisten quartz sand. After about 2 weeks, the seedlings reached a height of about 20 cm and were then transplanted into pots (10-cm diameter and 10-cm height, every 5 seedlings) containing 1/2 Hoagland nutrient solution. Seedlings were grown for 2 days to overcome any possible transplantation shock before starting the salt treatment. The salt treatment consisted of 300 mM NaCl in 1/2 Hoagland nutrient solution. Five seedlings of Jerusalem artichoke were distributed evenly in each pot. There were three replicate pots for each sampling time point.
In order to investigate the transcriptional responses in different tissues, we sampled roots and leaves during exposure to salt stress at 0 h, 6 h, 12 h, 24 h, and 48 h, with samples of roots named 0R, 6R, 12R, 24R, and 48R and samples of leaves named 0L, 6L, 12L, 24L, and 48L, respectively.

cDNA Library Establishment and Illumina RNA Sequencing
Total RNA of leaves and roots was extracted using a RNAiso Plus kit (Takara, Kyoto, Japan). The quality of RNA was detected by NanoDrop 2000 (Thermo, Waltham, MA, USA), and the integrity of the band was confirmed by gel electrophoresis (Bio-rad, Hercules, CA, USA). After digesting the DNA with DNase (Takara, Kyoto, Japan), we used magnetic beads with Oligo (dT) (Vdobiotech, Suzhou, China) to enrich eukaryotic mRNA. The RNA library was constructed using a TIANSeq Fast RNA Library Kit (Illumina) (Tiangen, Beijing, China). The library quality was assessed on an Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). Finally, samples were sequenced on a single lane of an Illumina HiSeq TM 2000 platform, and 125-bp paired-end reads were generated. To ensure the quality and reliability of data analysis, the raw data were filtered by removing reads with adapters, reads with unknown bases, and low-quality reads (Phred quality < 20). Software FastQC (Babraham Bioinformatics, Cambridge, UK) [33] was used to assess the quality of raw data before the data analysis. The low-quality reads were removed by a NGS QC Toolkit (version 2.3.3) [34]. The PCA used the default parameters of the R package FactoMineR [35,36].

De Novo Assembly
Because Jerusalem artichoke has no reference genome, de novo assembly of sequences was used for the follow-up analysis. The overlapping reads were connected into longer sequences; with their continuous extending, the reads were spliced into transcripts by the pair-end method using Trinity [37] (version 20131110). Subsequently, TGICL [38] was used for fast clustering and eliminating redundancy fragments to get final unigenes.

Unigene Annotation
Genetic similarity comparison was based on the BLAST method [39]. The databases RefSeq non-redundant protein sequences (NR) [40], SwissProt [41], Clusters of Orthologous Groups of Proteins (KOG) [42], Gene Ontology (GO) [43], and Kyoto Encyclopedia of Genes and Genomes (KEGG) [44] were used for comparison between unigenes and known genes, with the best similarity information chosen as the final annotation.

Expression Abundance and Selection of Differentially Expressed Genes (DEGs)
Gene expression abundance in each sample was estimated by bowtie2 [45] and eXpress [46]. The level of gene expression was calculated by the FPKM method [47] (Fragments Per kb per Million reads). Differential expression analysis was carried out using the DESeq [48] R package. The method of negative binomial distribution was used for identifying differentially expressed genes (DEGs) using the p-value < 0.05 level.

Time-Course Sequencing Data Analysis
MasSigPro R package was used for time-course sequencing data analysis (FPKM data) to identify significant differential gene expression profiles in RNA-seq [49]. The number of clusters was tested and found to have the optimal value of 9, and the other parameters used in the software were the default values.

Quantitative Real-Time q(RT)-PCR Analysis
Several DEGs that play important roles in salt stress were chosen to analyze for qRT-PCR analysis to check the consistency with the RNA-seq and qRT-PCR. Primer sequences for DEGs were designed by primer premier 6.0 (Table S1). The 25S rRNA was quantified as an internal control, and the 2 −∆∆Ct method [50] was used to analyze the differential expression. Total RNA was extracted by a RNAprep Pure Plant Kit (Polysaccharides-and Polyphenolics-rich) (Tiangen). A TaKaRa PrimeScript TM RT reagent Kit with gDNA Eraser (Perfect Real Time) was used for removing genomic DNA and reversing transcription. qRT-PCR was performed with an Applied Biosystems Step One Plus real-time PCR system with TaKaRa TB Green TM Premix Ex Taq TM (Tli RNaseH Plus). BIO-RAD Hard-Shell ® PCR Plates (96-well low-profile, semi-skirted) were used.

Phenotype Changes under High Salt Treatment
The phenotype of Jerusalem artichoke seedlings was recorded during the cultivation of plant material. Under the treatment with 300 mM NaCl, the change from 0 h to 12 h exposure was relatively small with roots darkening, stem beginning to bend, and leaves starting to wilt. After 24 h, the salt-induced changes were exacerbated. At 48 h, almost all leaves withered, and roots began to rot and become brittle ( Figure 1).

Evaluate and Annotation of All Unigenes
We evaluated the data after assembling the results of raw and clean reads, Q30 values, and GC content, indicating good quality of sequencing, with the data suitable for database construction, meeting the requirements for the subsequent analysis ( Figure S1B). The PCA showed large variability among the 10 samples ( Figure S1A). A total of 150,129 unigenes larger than 300 bp were used for statistical analysis. Most unigenes were in the 300-800 bp length range, with the number of unigenes decreasing with an increase in the length ( Figure S2A). The N50 value reached 1317, indicating the high quality of transcriptome assembly ( Figure S2B).
Due to a lack of Jerusalem artichoke genome information and other sequencing data, the function of a large number of unigenes is still unknown. Three databases were used for unigene annotations, including clusters of orthologous groups of proteins (KOG) (Figure 2A), Gene Ontology (GO) ( Figure 2B), and Kyoto Encyclopedia of Genes and Genomes (KEGG) ( Figure 2C). In KOG annotation, except for the "General function prediction only", the other most common one was "Signal transduction mechanisms", hinting at the importance of phytohormones in the response to salt stress ( Figure 2A). In the GO function database, the most common enrichment categories in "Biological Process" were "cellular process" and "metabolic process", in "Cellular Component" were "cell" and "cell part", and in "Molecular Function" were "binding" and "catalytic activity" ( Figure 2B). In the KEGG database, the groups with most terms were "Metabolism" and "Genetic Information Processing". The largest number of categories in "Genetic Information Processing" was in "Transcription" and "Translation", indicating that metabolic activity increased strongly in the high-salt treatment ( Figure 2C).

Differentially Expressed Genes
Except the two control (0 NaCl) groups of samples, there were eight experimental groups for identifying differentially expressed genes (DEGs). The genes with p-value ≤ 0.05 and absolute value of log 2 Fold change more than 2 were defined as differentially expressed genes (DEGs). The numbers of up-regulated and down-regulated unigenes between the experimental and control groups were identified ( Figure 3C). Due to the function of some DEGs being unknown, they needed to be filtered out; then, we chose the gene ontology (GO) classification. Finally, 9023 DEGs were identified and compared in Venn diagrams ( Figure 3A,B).

Time-Course Sequence Clustering and Expression Patterns of Genes Related to Plant Hormone Signal Transduction
Since the samples for RNA-seq were classified according to duration of salt treatment, we wanted to understand the similarities or differences in the expression patterns of the spliced unigenes. Finally, we divided the expression patterns of thousands of DEGs into nine categories for each roots and leaves using cluster analysis. After removing the lower expression branches, we finally obtained nine typical expression patterns ( Figure 4).
We chose the terms from the GO database for determining the similarity of the gene functions in the clusters. Interestingly, regardless of whether root or leaf samples were considered, redundant GO terms were positively correlated with the number of genes in the cluster, suggesting that functions were evenly distributed among modules, and modules did not have functional convergence ( Figure 5).  In order to further characterize the key genes in the clusters with the same expression pattern, we increased the screening threshold to p-value < 0.01 to observe the expression patterns and functional distribution (Figures 6 and 7).

Changes in Plant Hormone Signaling in the Transcriptome
Changes in plant hormones are critical in a variety of abiotic stress processes. In the KEGG enrichment analysis, we examined the enrichment of plant hormone signaling in each sample (Table 1). Compared with the no-stress controls (0R, 0L), the pathway of plant hormone signal transduction was enriched in almost all samples from NaCl-stressed seedlings. We screened all the unigenes that could be annotated to the pathways and classified them according to the annotations in the KOG database as the phytohormone type (auxin, cytokinin, gibberellin, abscisic acid, ethylene, brassinosteroid, jasmonic acid, and salicylic acid) and the participating pathways (Table 2) ( Figure 8A). Among the eight different plant hormone signal transduction pathways, unigenes were annotated to all types of phytohormones, indicating that the effects of salt stress on phytohormones were comprehensive and extensive. To explore further the expression patterns of these genes, we showed all the expression patterns of unigenes annotated to phytohormone signal transduction using heat maps ( Figure 8B).  To characterize the expression patterns of different phytohormones, we collected the expression data of all hormone-related genes and performed PCA analysis. The PCA results showed that the expression patterns of different hormones varied ( Figure 9A,B). In order to observe the relationship between different hormones, a co-occurrence network was constructed. Different hormones were in different positions in the network, and auxin, ethylene, and jasmonic acid were in the most critical positions in the network ( Figure 9C,D).

Quantitative Real-Time PCR Verification
In order to verify the entire transcriptome experiment and time-course analysis, we randomly selected 12 genes for verification by qRT-PCR. Similar to the sample set for transcriptome sequencing, we prepared samples of Jerusalem artichoke roots and leaves for gene expression at different durations of 300 mM NaCl stress treatment (0 h, 6 h, 12 h, 24 h, and 48 h). The results showed that the qRT-PCR experiment had good reproducibility, and the trends in gene expression were consistent with those revealed by the transcriptome sequencing ( Figures S3-S5), confirming the effectiveness of transcriptome sequencing ( Figure 10A). Every quantitative unigene had at least one functional coding sequence (CDS) ( Figure 10B).

Time-Course Expression Clustering and Functional Similarity
Traditional gene expression studies often focus on the difference between the control samples and the treatment samples. However, gene expression is strongly dependent on the duration of the treatment. Our time-course results showed that the expression of many genes was not linear with increased duration of salt stress (Figures 3, 6 and 7).
Many algorithms and programs have been developed to analyze time series of gene expression [51], enabling characterization of the dynamics of gene expression with increased treatment duration [49,52,53]. In the present study, we chose masigPro to cluster the expression data [49], and the parameter K = 9 to fully display the different expression patterns of all genes in either root or leaf samples. Among the nine expression clusters, some exhibited continuously rising or falling expression, whereas others reached the highest or lowest expression at a specific stress duration (6 h, 12 h or 24 h) ( Figures 5, 7 and 8).
The idea of predicting gene function through gene expression clustering to imply a possible relationship between temporal expression and gene function similarity [54] was not supported in the present study on Jerusalem artichoke ( Figure 5), but it is foreseeable that future studies of time-gradient transcriptomes may increase the number of sampling time points to explore changes in expression due to stress treatment in finer detail, even change in discrete time-point observations into continuous time observations.

Changes in Signal Transduction Pathways of Eight Plant Hormones
Plant hormones play a key role in the growth and development as well as in coping with various abiotic stresses [55][56][57][58][59]. In the present study, the pathways of plant hormone signal transduction under salt stress were enriched in eight comparison groups (based on plant organ and stress duration), showing complexities of the phytohormonal responses to salt stress (Figure 3) [60]. In this research, genes related to signal transduction of all eight types of plant hormones were involved, suggesting that phytohormones regulate the transcription dynamics in Jerusalem artichoke in the early stage of high salt stress (Figure 11). Among the plant hormones, ABA is one of the most widely studied ones. PYL (from the abscisic acid receptor PYR/PYL family) is a protein that mediates ABA signaling, acting as a receptor for ABA that regulates a series of downstream responses. In the present study, seven of the nine genes belonging to the PYL family were highly expressed in roots at 0 h, and the expression was decreased with increasing duration of salt stress, indicating that the PYL expression was suppressed by the salt treatment. PP2C (Type 2C protein phosphatases) is also a key link in the ABA signaling pathway, with PYR/PYL controlling ABA signaling by inhibiting PP2C [61]. In our results, it was found that the expression of PYL in the early stage of root stress (0-6 h) was opposite to the expression of PP2C, which also confirmed the existence of this inhibitory effect, which gradually disappeared over time. Hence, the inhibition of PP2C occurred in the early stage of salt stress ( Figure 8B), which is consistent with the results of the Arabidopsis salt-stress expression profile [32].
In the auxin signaling pathway, GH3 (Gretchen Hagen3) was highly expressed only in roots in the early stage of salt stress (6 h) ( Figure 8B). The GH3 gene is involved in auxin homeostasis by catalyzing auxin conjugates and binding free indole-3-acetic acid (IAA) to amino acids. The GH3 mRNA has been shown to be one of the most relevant early molecular markers of primary auxin response [62,63]. The GH3 gene family in maize responded to several abiotic stresses and stress-related phytohormones [64].
Ethylene has long been recognized as a stress-related phytohormone that regulates different levels of salt response, including membrane receptors, components in the cytoplasm, and nuclear transcription factors [65], with the ethylene receptor (ETR) being the first component of ethylene signaling. Overexpression of the Nicotiana tabacum type II ETR homologue NTHK1 in Arabidopsis increased salt sensitivity, and the mutants of EIN2 (ethylene-insensitive protein 2) also had increased sensitivity to salt stress [66], suggesting that EIN and ETR were important in a response to salt stress in the ethylene signaling pathway. In non-model plant mulberry (Morus sp.), the expression of the gene encoding EIN3 protein was induced by salt stress. Over-expression of that gene in Arabidopsis enhanced its salt resistance [67]. Although there are many reports on ethylene and salt stress, the temporal and tissue-specific expression of genes involved in ethylene signal transduction under salt stress is still poorly understood. In the present study, the expression of most genes in the ethylene signaling pathway was increased in roots after 24 h of salt exposure. In leaves, except for ERF, most genes involved in ethylene signaling were expressed in the early stage of stress treatment (6 h), and there was a tendency of increased expression with the duration of salt stress (peaking for most genes at 24 h) ( Figure 8B).
The expression of Pohlia nutans L. jasmonate ZIM-domain gene PnJAZ1 was induced by salt stress [68]. Transgenic Arabidopsis thaliana overexpressing PnJAZ1 had increased tolerance to salt stress and decreased ABA sensitivity during early development [69]. In the present study, up to 16 JAZ genes were annotated into the jasmonic acid signaling pathway ( Figure 8B). The Arabidopsis JIN1/MYC2 mutant jin1 (jasmonate insensitive 1) had significantly lower resistance to salt stress and lower proline and sugar content than the wild type. The concentration of flavonoids demonstrated that MYC2 was involved in controlling the plant protection systems under salt stress [70]. In the present study, the expression patterns of differentially expressed genes (DEG) in the jasmonic acid signaling pathway were similar, being highly expressed in roots in the early stage (6 h) of salt stress but not at other times ( Figure 8B).
The expression of genes related to other four plant hormones' (brassinosteroids, gibberellins, cytokinins, and salicylic acid) signal transduction did not show distinct patterns ( Figure 8B).

Time-Course Expression of Phytohormone Signal Transduction Genes in Jerusalem Artichoke
Phytohormonal changes were large in the initial stage of plant stress. The results of the enrichment analysis (Table 2) were important in explaining the role of phytohormone signal transduction during exposure to high salt stress. We characterized the expression patterns of the genes related to phytohormonal signal transduction in the time course of salt-stress exposure. This type of study is rare [71,72].
In our salt-stress study of Jerusalem artichoke, the expression pattern of the phytohormone signal transduction-related genes was highly dependent on the type of phytohormone. Similar findings were reported in Arabidopsis subjected to moderate dehydration stress [73]. The complex signal transduction network generated by the interaction of different phy-tohormones was also found under copper toxicity treatment, explaining the differential expression patterns of various phytohormones [74]. The variable temporal expression patterns of the phytohormone metabolism genes were reported during the process of water absorption in Arabidopsis seeds [75]. The above-mentioned literature, together with our own findings, indicate the relevance of characterizing the time-course of the expression of phytohormone-related genes and functions.

Conclusions
This was the first time-series analysis of the Jerusalem artichoke transcriptome under salt treatments. The upper limit of salt tolerance of Helianthus tuberosus in laboratory culture is 300 mM. Under the impact of a high-salinity (300 mM) environment, Jerusalem artichoke in the seedling stage struggled to survive for a long time, and the phenotype was severe in the short term. In the first 48 h of salt-stress response, Jerusalem artichoke transcription changes drastically, and the subsequent transcription changes will tend to be stable. Discovering the early transcription changes is critical to understanding the salt tolerance mechanism. Based on the expression of genes on the time scale, we found that the distribution of gene functions in time is relatively even. Upregulation of the phytohormone signal transduction had a crucial role in the response of Jerusalem artichoke seedlings to salt stress, with genes related to each phytohormone showing a similar expression pattern during the first 48 h of exposure to salt stress. The genes in four of the eight hormone (abscisic acid, auxin, ethylene, and jasmonic acid) pathways had the most obvious change pattern. The revealed dynamics of genes related to plant hormone signal transduction under salt stress represent a basis for further characterization of the downstream gene regulation.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/agronomy12010163/s1, Table S1. Primers used in qRT-PCR; Figure S1. PCA (Principal Component Analysis) of 5 root (0R, 6R, 12R, 24R, 48R) and 5 leaf (0L, 6L, 12L, 24L, 48L) samples (A). Details of the raw reads, clean reads, Q30 and GC content (B); Figure S2. Length distribution of unigenes meeting conditions (A); Information on assembled unigenes. N50: the sequence length of the shortest scaffold at 50% of the total genome length, representing the quality of assembly (B); Figure S3. The qRT-PCR validation of 4 unigenes (CL102748Conitg1, CL82936Contig1, CL862Contig2, CL42555Contig2). The histograms in each plot represent the results of real-time qRT-PCR, and the line graph represents the results of RNA-seq. The test method for the difference of qRT-PCR results at different times used one-way ANOVA (p-value < 0.05). Differences among samples are represented by different letters; Figure S4. The qRT-PCR validation of 4 unigenes (Comp37929_c0_seq2_F, CL961Contig1, CL70841Contig1, CL133346Contig1). The histograms in each plot represent the results of real-time qRT-PCR, and the line graph represents the results of RNA-seq. The test method for the difference of qRT-PCR results at different times used one-way ANOVA (p-value < 0.05). Differences among samples are represented by different letters; Figure S5. The qRT-PCR validation of 4 unigenes (CL5042Contig1, CL69482Contig1, CL29820Contig1, CL5184Contig2). The histograms in each plot represent the results of real-time qRT-PCR, and the line graph represents the results of RNA-seq. The test method for the difference of qRT-PCR results at different times used one-way ANOVA (p-value < 0.05). Differences among samples are represented by different letters.
Author Contributions: Z.Z., X.L. and Y.Y. designed the experiments. Y.Y., X.G., W.R. and Y.Y. performed the experiments. Y.Y., J.W., Z.Z. and X.L. analyzed the data. Y.Y., Z.Z., X.L. and Z.R. wrote and revised the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The raw data of RNA-seq have been uploaded to NCBI. BioProject ID is PRJNA793515.