Genome-Wide Identification of GRAS Gene Family and Their Responses to Abiotic Stress in Medicago sativa

Alfalfa (Medicago sativa) is a high-quality legume forage crop worldwide, and alfalfa production is often threatened by abiotic environmental stresses. GRAS proteins are important transcription factors that play a vital role in plant development, as well as in response to environmental stress. In this study, the availability of alfalfa genome “Zhongmu No.1” allowed us to identify 51 GRAS family members, i.e., MsGRAS. MsGRAS proteins could be classified into nine subgroups with distinct conserved domains, and tandem and segmental duplications were observed as an expansion strategy of this gene family. In RNA-Seq analysis, 14 MsGRAS genes were not expressed in the leaf or root, 6 GRAS genes in 3 differentially expressed gene clusters were involved in the salinity stress response in the leaf. Moreover, qRT-PCR results confirmed that MsGRAS51 expression was induced under drought stress and hormone treatments (ABA, GA and IAA) but down-regulated in salinity stress. Collectively, our genome-wide characterization, evolutionary, and expression analysis suggested that the MsGRAS proteins might play crucial roles in response to abiotic stresses and hormonal cues in alfalfa. For the breeding of alfalfa, it provided important information on stress resistance and functional studies on MsGRAS and hormone signaling.


Introduction
Transcription factors (TFs) as regulatory proteins could bind to specific DNA sequences (cis-acting elements), which might be located in the promoters of target genes, and play a vital role in plant development, as well as in response to abiotic stresses, such as drought, salt, chilling and heat [1]. The GRAS transcription factors are proposed to be plant-specific regulation proteins, which are supported by their possession of certain structural features [2]. The GRAS family was defined, after GAI (gibberellic acid insensitive), RGA (repressor of GA1-3 mutant), and scarecrow (SCR) were identified as family members. GAI and RGA play a vital role in the gibberellin-dependent signal transduction process, and SCR is involved in regulating the radial tissue differentiation in roots [3][4][5]. Typically, the GRAS protein sequences consist of 400-700 amino acids (aa), including a highly conserved C-terminal and a variable N-terminal. Its C-terminal is composed of five conserved domains, namely SAW, LRI, LRII, PFYRE, and VHIID [2,6]. The N-terminal is the functional component of the GRAS protein, as its IDRs (intrinsically disordered regions) could bind with different other proteins [7].
We carried out a collinearity analysis on MsGRAS genes with those in a monocotyledonous plant of rice (Oryza sativa) and three dicotyledonous plants of Arabidopsis, soybean, and M. truncatula (Figure 2), to further understand the collinearity of the MsGRAS genes. MsGRAS genes displayed various numbers of syntenic lines with four species. A total of 36 MsGRAS genes showed a syntenic relationship with those in soybean, and 28 corresponding orthologs were identified between M. sativa and M. truncatula. Meanwhile, only 6 and 12 MsGRAS genes showed syntenic relationships with those in O. sativa and A. thaliana, respectively ( Figure 2 and Table S2).
Sequence analyses revealed that the size of the MsGRAS proteins vary greatly in length, from 98 aa (MsGRAS24) to 1002 aa (MsGRAS12). The proteins molecular weight ranged from 23.44 kD (MsGRAS49) to 112.7 kD (MsGRAS12), and the isoelectric point (pI) ranged from 4.77 (MsGRAS45) to 9.24 (MsGRAS15). Subcellular localization prediction by WoLF PSORT tool showed that 33 MsGRAS proteins are located in the nucleus, 12 in the cytoplasmic, and 6 in the chloroplast (Table 1), while TargetP prediction results were all "Other", rather than chloroplast or mitochondria. The results showed that MsGRAS protein might function in diverse microcellular environments.   The online MEME search tool was used to further look into the conversation of MsGRAS proteins, and 10 distinct conserved motifs (named as motif 1-10) were found (Figure 3a,b and Table S3). Based on the alignment of predicted GRAS domain sequences, we found that most of the MsGRAS proteins (94%) contained motifs of 4, 6, and 7, and 2 copies of motif 7 were separated by motifs of 2 and 3 in MsGRAS18. Four protein pairs of MsGRAS2/3, MsGRAS7/8, MsGRAS10/11, and MsGRAS39/40 shared identical motif patterns, and MsGRAS24 held only three motifs of 1, 5, and 10. Multiple sequence alignments showed that most of the MsGRAS proteins possess the N-terminal and the highly conserved C-terminal, and the conserved domains located in the C-terminal included VHIID, LHRI, and SAW. Among the whole MsGRAS proteins, 47 MsGRAS contain VHIID and LHRI domains, and 49 proteins contain SAW domains ( Figure 4 and Figures S1-S3). The motifs were located in the four GRAS domains, including LHRI (motif1, motif5, motif6, motif10), VHIID (motif3, motif7), PFYRE (motif2, motif8), and SAW (motif4, motif9), and shared across almost all MsGRAS members. We further found that MsGRAS24 with the gene ID of MsG0480021771.01.T01 is wrongly annotated in the alfalfa genome reference of "Zhongmu No.1". Based on our unpublished PacBio sequencing dataset for a pool of root, stem, leaf, etc., in alfalfa "Zhongmu No.1", one full-length transcript was retrieved for this gene (Table S1), and based on the genomic coordinates of this transcript, its gene annotation and structure were updated without long introns ( Figure 3c).

Transcriptome Analysis of the Alfalfa Leaves under Abiotic Stress
To investigate the expression profiles of MsGRAS members under abiotic stresses (salt and drought), we downloaded publicly available leaf and root transcriptome data from NCBI (Bioproject PRJNA657410 and PRJNA525327) [35] and conducted RNA-seq analysis. The alfalfa plant was treated with salt stress for 3 (L3 for leaf and R3 for root) and 27 h (L27 for leaf and R27 for root) with the control (L0 for leaf and R0 for root), while drought stress was also applied for roots under drought stress (RD). In root tissue, 35 MsGRAS genes were expressed in R0 (Table S4), and the expression of 4 and 5 MsGRAS genes changed significantly in R3 and R27, respectively, including MsGRAS51/13. In leaf tissue, 29 MsGRAS genes were expressed in L0 (Table S4), and the expression levels of 17 and 19 MsGRAS genes increased in L3 and L27, respectively ( Figure 6, Table S4). Actually, 14 MsGRAS genes were not expressed in both root and leaf. Under the drought stress, 13 MsGRAS genes were insensitive, and about 41% of MsGRAS genes showed significant differences, of which eight genes belonged to the PAT1 subgroup. Furthermore, we performed a differential expression analysis of all genes in alfalfa leaves under salinity stress compared with the control groups. A set of differentially expressed genes (DEGs) for each salinity stress treatment were determined according to their significance with a cutoff of p-value < 0.05 and |log2(fold change)| > 1 (Table S5 and Figure 7a). The result showed that a total of 1090, 1975, and 929 DEGs, in response to NaCl treatments, were detected in the three comparisons (L0 vs. L3, L0 vs. L27, L3 vs. L27), respectively. The DEGs with down-regulated ones in L0 vs. L27 were significantly more than those in L0 vs. L3. Compared to the control, there were 526 up-regulated and 564 down-regulated DEGs under 3 h salinity stress, whereas 837 up-regulated and 1138 downregulated DEGs in response to 27 h NaCl treatment in alfalfa leaves. A total of 929 DEGs were detected in L3 vs. L27, including 534 up-regulated and 395 down-regulated ones.
The Venn diagram further showed that 141 DEGs were common in the three comparisons ( Figure 7b). These results suggested that the alfalfa plant was more sensitive to salinity stress for a longer time, e.g., 27 h in this case, with more DEGs. To trace the partners of MsGRAS genes in salinity response, we analyzed the coexpression of these DEGs and divided them into 12 clusters. A total six MsGRAS genes were placed in three clusters, i.e., MsGRAS37 in cluster 1; MsGRAS5, MsGRAS13, and MsGRAS51 in cluster 5; MsGRAS27 and MsGRAS50 in cluster 7 (Figure 7c). In cluster 1, the gene expression decreased rapidly at 3 h of salt stress treatment and maintained at 27 h. The gene expression in cluster 5 did not change significantly at 3 h and showed an obvious down-regulation trend at 27 h, while the gene expression pattern in cluster 7 is up at 27 h.
To identify functions of salinity-responsive partners of MsGRAS genes in alfalfa, we conducted GO term enrichment analysis on DEGs in cluster 1, cluster 5, and cluster 7 in leaf. In the GO enrichment graph, the p-value is sorted from small to large to indicate the significance of the enrichment degree (Figure 8a). The GO enrichment results showed that the majority of DEGs in cluster 1 were involved in molecular function categories, including 'catalytic activity', 'transferase activity', 'transferase activity transferring glycosyl groups', and 'DNA binding'. The DEGs in cluster 5 were involved mostly in molecular functions, including 'catalytic activity', 'oxidoreductase activity', 'cofactor binding', and 'transferase activity'. Those in cluster 7 were mainly enriched in terms of 'lysine methyltransferase activity', 'transcription factor activity', and 'protein heterodimerization activity'. The expression of MsGRAS5, MsGRAS13, MsGRAS37, and MsGRAS51, which belong to the DELLA and PAT1 subfamilies, decreased significantly at 27 h (Figure 8d), which indicated the role of DELLA and PAT1 proteins in response to salt stress.

Validation of Expression Profiles of MsGRAS Genes under Abiotic Stress and Hormone Treatments
Among all 51 MsGRAS genes, expressions of 6 (MsGRAS5, MsGRAS13, MsGRAS27, MsGRAS37, MsGRAS50, and MsGRAS51) were highly induced by salinity and drought stresses. We designed primers for the six selected genes and conducted qRT-PCR experiments in leaf tissues (Table S6). The gene expression profiles of the six selected genes were generally consistent with the RNA-seq data (Figure 8d and Table S7). The results showed that upon salinity treatment, MsGRAS50 expression was significantly induced (p < 0.01), while the other four MsGRAS genes' expression was down-regulated significantly. Under drought stress, the expression of MsGRAS51 increased gradually. After 72 h of treatment, MsGRAS51 expression was about 20 times that of the control, implying its putative roles in stress tolerance. Interestingly, the expression of three genes (MsGRAS5, MsGRAS13, and MsGRAS50) were only up-regulated at 12 h. However, MsGRAS27 showed no significant influence in both salt and drought treatment (Figure 9).
Previous studies have demonstrated that GRAS genes are key regulators of germination [21], seedling growth [2], lateral root development, and shoot meristem maintenance [6]. Plant hormones have been widely reported for their roles in the regulation of various aspects of plant growth. Upon drought and salinity stress, gene expression profiles showed that MsGRAS genes play a vital role in response to abiotic stress. We further verified the expression profiles of 6 selected MsGRAS genes in alfalfa seedlings under exogenous ABA, GA, and IAA treatments by qRT-PCR. We found that the hormone treatments resulted in a wide variety of MsGRAS gene expression profiles ( Figure 10, Table S8). It showed that MsGRAS5 was not responsive to ABA, and the other five genes' expression, with one up-and four down-regulated, was significantly different from the control. Four MsGRAS genes showed significantly different expressions after treatments of both GA and IAA. The up-regulated expression levels of MsGRAS51 responded to all three hormones. It suggested that the MsGRAS TFs played vital roles in response to ABA, GA, and IAA stress signaling. Under GA treatment for 72 h, MsGRAS50 was significantly up-regulated, with >800 folds of the control. MsGRAS5 as DELLA protein showed a significant response to gibberellin treatment and no expression at all at 24 h. The expression levels of six MsGRAS genes varied obviously in response to hormone treatments, suggesting that the MsGRAS genes have been involved in different pathways in alfalfa seedlings.

Duplication of MsGRAS Genes
In this investigation, the identified MsGRAS genes are unevenly distributed on eight chromosomes of the alfalfa genome. Notably, Chr2 and Chr4 contained most of MsGRAS genes (20 genes), which mainly originated from all of the nine subfamilies, and Chr1 is the "cold" region, including only one MsGRAS gene (Figure 1a). In M. truncatula, Chr2 and Chr4 are also the "hot" regions, with 15 MtGRAS genes [13]. Furthermore, based on the distribution of MsGRAS in chromosomes, we discovered the gene duplication events. Gene loss and replication events greatly promote the evolution of species, providing raw materials for generating novel gene functions and expanding gene families [43]. Most plants have experienced one or more rounds of ancient polyploidy, such as Arabidopsis, which had undergone two recent whole-genome duplications [44]. It is speculated that gene replication may be one of the main factors for the number variation of GRAS family members. We discovered 12 pairs of MsGRAS genes as evidence of tandem duplication events on five chromosomes (Chr2, Chr4, Chr5, Chr7, and Chr8), among which six belonged to the LISCL subfamily, and three segmental duplication events for ten MsGRAS genes (Figure 1b). It seems that tandem duplication contributed more to M. sativa GRAS expansion than segmental duplication. For example, LISCL subfamily expansion in alfalfa is also found in M. truncatula [13]. Segmental duplication and tandem duplication are considered to represent principal evolutionary patterns in plants. Segmental duplications are produced by polyploidy and chromosomal rearrangement because of numerous duplicated chromosomal blocks in plant genomes [45]. Interestingly, almost all genes in MsGRAS families have no introns. In the GRAS gene family, the higher percentage of intronless genes, the closer the evolutionary relationship of the GRAS members [46]. Intron deficiency is common in other large gene families, for example, DEAD-box RNA helicases, small auxin-up RNAs (SAUR) gene family, and the F-box transcription factor family [47][48][49]. Intronless genes are also common in prokaryote genomes, but there are three explanations for the appearance of intronless genes in eukaryote genomes, i.e., horizontal gene transfer, intronless gene reverse transcription, and intronless gene replication [46]. In addition, the GRAS domain proteins were discovered in bacteria [50], and thus, horizontal gene transfer from ancient prokaryote genomes might be a possibility to explain the abundant intronless genes in the plant GRAS gene family.

Expression and Function of MsGRAS Genes
The GRAS members in the same subfamily commonly exhibit functional similarities. Structural analysis provides a clue to locate which subgroup of GRAS is of ancient origin [51]. All of the ten predicted alfalfa DELLA proteins (MsGRAS5, MsGRAS13, Ms-GRAS15, MsGRAS26, MsGRAS29, MsGRAS30, MsGRAS31, MsGRAS32, MsGRAS33, and MsGRAS38) were classified into the DELLA subgroup, together with the Arabidopsis DELLA proteins of AtGAI, AtRGA, AtRGL1, AtRGL2, and AtRGL3. In Arabidopsis, the DELLA family consists of GAI, RGA, RGL1, RGL2, and RGL3 [4,21,23,52]. To date, the functions of the AtGAI and AtRGA proteins have been clearly illuminated in Arabidopsis, and DELLA protein, as a putative transcriptional regulator, is a subbranch of the GRAS gene family, playing a negative regulatory role in gibberellin (GA) signal transduction [53]. DELLA proteins restrain the development of plants, whereas GA relieves plants of DELLA-mediated growth restraint [54,55]. After GA treatment, the expression of MsGRAS5 significantly changed, especially without expression at 24 h, indicating its potential role in gibberellin regulation in alfalfa.
It is known that GRAS proteins might interact with each other in a pathway. The GRAS family transcription factors of SCR and SHR are involved in root tissue formation, and the two subfamilies AtSHR and AtSCR, were considered to function importantly in maintaining stem and root cell meristem [12,56]. As a positive regulator, SCL3 acts to integrate and maintain a functional GA pathway by eliminating the negative regulation of DELLA to control the dynamic balance of GA in root development for cell elongation [57]. In this study, the SCL3 subfamily contained four MsGRAS proteins (MsGRAS12, 24, 27, and 45), which is implying a similar function in root development. In addition, LISCLs are of the largest GRAS subfamily in alfalfa and involved in the regulation of microsporogenesis in Lilium longiflorum [43], while LAS proteins act as branching regulators, which is required to establish axillary meristems at a distance to the SAM during vegetative development [24,58,59].

Possible Functions of MsGRAS Members
The gene expression profiles could provide key clues to predict their possible functions. As a broad-spectrum phytohormone, ABA is involved not only in regulating plant growth and development but also in coordinating in many aspects of stress signal transduction pathways during abiotic stresses. ABA mediates many stress responses and activation of genes involved in tolerance adjustment [60]. When plants are subjected to biotic and abiotic stress, cells are prone to oxidative stress, which is due to the increase of the content of reactive oxygen species (ROS) [61]. The phytohormone GA is known to play an important role in regulating a diverse array of developmental processes, such as seed development and germination, organ elongation, and control of flowering time. The GRAS subfamily has been reported to function as repressors of gibberellin (GA) signaling [2]. A previous study reported that DELLA proteins thus repressed ROS accumulation, interrupts ROS-induced cell death, and hence enhance plant biotic and abiotic stress tolerance [62].
Recent studies have shown that GRAS genes have multiple functions in plant abiotic stress tolerance. The rice GRAS transcription factor gene OsGRAS23 was induced by osmotic stress, and overexpression of this gene enhanced the drought resistance in transgenic rice plants by regulating the expression of stress-responsive genes [63]. In M. sativa, MsGRAS37 and MsGRAS51 showed a high expression level under PEG treatment. Moreover, the expression levels of several MsGRAS genes (e.g., MsGRAS5/13/50/51) were dramatically regulated by salinity stresses. It suggested they played critical roles in the abiotic stress tolerance in M. sativa. In our study, the expression of MsGRAS51 significantly responded to all the treatments of hormone and stress, indicating a potential regulatory role in alfalfa development and stress resistance.

Plant Materials
In this study, the seeds of M. For the expression analyses of GRAS genes' response to different treatments, the seven-day-old seedlings were transferred into the flasks with 1/2 MS liquid medium and grew in a controlled growth chamber under 16/8 h light/dark regime at 20 • C. After the third leaf was fully expanded, the plants were supplied with 1 mmol·L −1 IAA, 1 mmol·L −1 GA3, 1 mmol·L −1 ABA, 15% PEG6000, and 300 mmol·L −1 NaCl to the 1/2 MS liquid medium, respectively. For each treatment, the whole plant was collected at 0, 2, 24, and 72 h with a triplicate.

RNA Extraction and Gene Expression Analysis
The total plant RNAs were extracted using the HUA YUE YANG total RNA Extraction kit (Beijing, China), according to the manufacturer's instructions. RNA was used as a template for the synthesis of first-strand cDNA. Complementary DNA was generated with reverse transcriptase (TransGen Biotech, Beijing, China) at 42 • C for 15 min, 85 • C for 5 s. Quantitative real-time PCR (qRT-PCR) was conducted according to the instructions of 2 × RealStar Green Fast Mixture (GeneStar, Beijing, China) on CFX96 Touch TM RT-PCR system (Biorad, Los Angeles, CA, USA). The gene-specific primer sequences for qRT-PCR determination are provided in Table S6. Ms-Actin was used as the internal control. Three technical repetitions for each sample were performed, and the relative expression data were calculated according to the 2 −∆∆CT method. Student t-test was performed by Graphpad Prism 7 (https://www.graphpad.com/scientific software/prism/ accessed on 5 April 2021). All the error bars were standard deviations (SD) from three biological replicates.

Identification of GRAS Genes in M. sativa
The M. sativa "Zhongmu No.1" genome information [31] was downloaded from the figshare website (https://figshare.com/articles/dataset/Medicago_sativa_genome_ and_annotation_files/ accessed on 11 September 2020). The 32 GRAS protein sequences of Arabidopsis retrieved from the TAIR (https://www.arabidopsis.org/ accessed on 11 September 2020) were used as a query to obtain the possible GRAS proteins in M. sativa genome by BlastP search with a cutoff E-value of 1.0 × 10 −10 . Furthermore, the GRAS protein sequences of M. sativa were aligned using the HMM model of the HMMER3.0 (http: //hmmer.org/ accessed on 11 September 2020) (National Institutes of Health, Bethesda, MD, USA), with an E-value of 0.001 and the removal of redundant sequences. The Hidden Markov Model (HMM) profile (PF03514.13) was used to identify the putative GRAS domain in the Pfam database (http://pfam.xfam.org/ accessed on 11 September 2020) [64].

Chromosomal Distribution and Gene Duplication Events Analysis
The chromosomal locations of 51 GRAS genes were obtained from the M. sativa genomic annotation file GFF3 (general feature format). We draw the chromosomal distribution image of MsGRAS genes through the TBtools software [65]. The 51 putative GRAS genes were renamed according to their chromosomal locations. The detection and identification of gene duplication events in MsGRAS genes were performed using multiple collinear scanning toolkits (MCScanX) [66] with an E-value set to 10 −5 .

Sequence Analyses and Structural Characterization of MsGRAS Genes
The tools from the ExPASy website (https://web.expasy.org/protparam/ accessed on 11 September 2020) (SIB Swiss Institute of Bioinformatics) were used to analyze the basic physical and chemical characteristics, such as sequence length, molecular weight, and isoelectric point of the identified GRAS proteins. The subcellular localization of MsGRAS was predicted on the WoLF PSORT (http://wolfpsort.org accessed on 10 July 2021). The exon-intron structures of MsGRAS genes were drawn based on the CDS and the corresponding full-length sequence by using TBtools software [65]. The conserved motifs in MsGRAS proteins were analyzed through the MEME software [67]. The parameters were set as follows: motif width as 15-50 amino acids (aa) and the number of motifs as 10.

Phylogenetic and Collinearity Analysis of MsGRAS Genes
We downloaded the GRAS proteins of Arabidopsis and M. truncatula from the phytozome (https://phytozome.jgi.doe.gov/pz/portal.html accessed on 5 October 2020), which were used for the phylogenetic analysis of GRAS proteins in plants. Multiple sequence alignments of the 142 GRAS genes were performed using the ClustalX software (https://clustalx.software.informer.com/2.1/# version 2.1, accessed on 5 October 2020) with default parameters. The phylogenetic tree based on the alignment was constructed using MEGA-X [34] by neighbor-joining (NJ) method, and the reliability of the obtained trees was assessed with a bootstrap value of 1000. The phylogenetic tree was visualized and illustrated using the online tool EvolView [68] (http://www.evolgenius.info/evolview/# version 3.0, accessed on 5 October 2020). The syntenic relationship between MsGRAS genes and GRAS genes from Arabidopsis, rice, M. truncatula, and soybean were determined by using Dual Synteny Plotter (Plant Genome Mapping Laboratory, University of Georgia, Athens, GA, USA, http://chibba.pgml.uga.edu/mcscan2/# version 2.0, accessed on 5 October 2020) software and was visualized by TBtools software (HuaZhong Agricultural University, Wuhan, China).

Conclusions
We screened 51 MsGRAS genes in the alfalfa genome of "Zhongmu No.1" and grouped them into nine subfamilies, which all shared the conversed GRAS domain. We also identified the conserved domain patterns and phylogenetic relationships of these MsGRAS proteins. Our current systematic studies provided gene expression atlas for MsGRAS genes. These results provide comprehensive information towards better understanding the molecular basis and impact of GRAS gene families on the response of alfalfa to treatments of hormones and stress. Considering that the detailed regulatory mechanisms underlining GRAS genes are still not clear in M. sativa, it will be of great interest to elucidate individual MsGRAS genes in the near future. For instance, MsGRAS51 might contribute to abiotic stress tolerance in M. sativa, and MsGRAS50 may play an important role in the gibberellin signal transduction pathway. Nevertheless, it provided important information about the GRAS family and a framework for breeding on stress-resistance and hormone signaling transduction in alfalfa.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/ijms22147729/s1, Figure S1: LHRI domain of 51 MsGRAS proteins based on sequence alignment. The most conserved GRAS domain of LHRI was boxed in red; Figure S2: SAW domain of 51 MsGRAS proteins based on sequence alignment. The most conserved GRAS domain of SAW was boxed in red; Figure S3: PFYRE domain of 51 MsGRAS proteins based on sequence alignment. The most conserved GRAS domain of PFYRE was boxed in red; Table S1: Nucleotide and amino acid sequences of MsGRAS; Table S2: The collinearity relationships of the MsGRAS genes in four species; Table S3: Analysis and distribution of conserved motifs in MsGRAS proteins; Table S4: The relative expression level of 51 MsGRAS genes under various treatments of abiotic stress and treatment time; Table S5: Gene expression level of 51 MsGRAS genes in root and leaf tissues under treatments of salt stress; Table S6: The six gene-specific primer pairs used for qRT-PCR validation; Table S7: Expression levels of six MsGRAS genes under drought, salt, exogenous ABA, GA and IAA treatments by qRT-PCR; Table S8: Expression data of 51 MsGRAS genes in public RNA-seq datasets under stress treatments.
Author Contributions: H.Z. performed stress treatments, data analysis, gene expression analysis, and wrote the manuscript; X.L. prepared the figures. X.W., M.S. and R.S. collected the data; S.J. and P.M. designed the experiments, supervised this work, and edited the manuscript. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: RNA-seq raw data used in this study was downloaded from the SRA database in NCBI with the Bioproject accession numbers of PRJNA657410 and PRJNA525327.