Leaf transcriptome analysis of Medicago ruthenica , revealing its response and adaptive strategy to drought stress and rehydration

Background: Recently, drought stress has brought tremendous loss on the production of agriculture and animal husbandry. In realistic production, plants are often in cyclic wet-dry environment. Therefore, the factors that affect the final yield of plants in adversity including the resistance and tolerance to drought and the ability of plants to resume from the previous damage after rehydration. So it’s necessary for us to study the response and adaptive strategies of plants to drought and rehydration. Generally, the yield of herbage with strong resistance is relatively low. However, Medicago ruthenica (L.)cv.Zhilixing has the advantages of strong resistance and high yield concurrently. This made it can be used for raising livestock, natural grassland improvement, as a good parent for breeding and a new and high quality resource of stress resistance genes. Now, there are still many problems need to be solved when compared with other important legume forages. Therefore, we analyzed the changes of Medicago ruthenica (L.)cv.Zhilixing on transcription level under drought stress and rehydration, explored its phased response strategies. Results: We obtained 191 DEGs in drought stress, and the three treatments has 43 DEGs in common. Galactose metabolism, Starch and sucrose metabolism, Arginine and proline metabolism, TCA cycle, Photosynthesis-antenna proteins, were involved in the adaptation of Medicago ruthenica to 9 days of drought stress. The regulation of Arginine and proline metabolism, Cysteine and methionine metabolism, Photosynthesis-antenna proteins, Ascorbate and aldarate metabolism were conducive to the resistance of Medicago ruthenica to severe drought stress. The regulation of Starch and sucrose metabolism, Flavonoid biosynthesis, Valine, leucine and isoleucine degradation, Circadian rhythm-plant was beneficial to the post drought recovery of Medicago ruthenica . Conclusions: We preliminarily analyzed the adaptation mechanism of the plant under different drought and rehydration conditions. Medicago ruthenica (L.)cv.Zhilixing adopts different strategies to adapt to different degrees of drought stress and rehydration. The research discovered the genes that can be used as candidate genes to improve stress resistance and drought adaptability of plants. Our transcriptome data dramatically enriches the resources of stress resistance genes. It can provide theoretical support for further adaptation mechanism research of the plant under different drought and rehydration conditions. data and basis for further revealing the molecular mechanism of Medicago ruthenicus adaptation to drought and rehydration.

found that the leaves were elliptic with serrated edges, and the upper epidermis was sparsely pubescent, while the lower epidermis was more abundant than the upper epidermis. The stomatal depression was deeper than that of alfalfa, and the stomatal density, the compactness of palisade tissue and spongy tissue, and the layer number of palisade tissue were also higher than that of alfalfa, which was conducive to photosynthesis and prevented water evapotranspiration in drought.
Meanwhile, compared with alfalfa, the leaf veins were more obvious and dense, the conducting tissue was developed, and the intercellular space was smaller. Morphological and anatomical analysis showed that the erect type of Medicago ruthenica was more drought resistant than alfalfa.
Furthermore, this variety is not only a good parent for breeding, but also a good gene resource for improving the stress resistance of alfalfa and other forages. It is better than alfalfa in adaptability and resistance. Therefore, it is of great significance on the research, development and utilization of the plant.
Using transcriptome sequencing technology to study the expression of the whole genome in plants under drought stress is helpful to reveal the molecular regulation mechanism of plant response to drought. It is of vital importance to identify the key genes in various metabolic pathways, to explore new genes that related to drought resistance and to construct transcriptional regulatory networks in plant response to stress [12]. RNA-seq has been used to study the response of Medicago ruthenica to single stress for 3 hours (ABA, low temperature, salt stress, etc.) [13].
However, there is no report on the transcriptional changes of the forage under different levels of drought stress and rehydration. Therefore, the experiment aims to study the variations of the molecular level of Medicago ruthenica (L.) cv. Zhilixing under different intensity of drought stress and rehydration. It's benefit to understand the molecular response of the plant to drought stress and rehydration, and also to provide a reference for the later study of the phased response strategy of Medicago ruthenica to drought stress and rehydration. The research will provide a new idea for improving the drought resistance and drought adaptability of alfalfa and other forages. And also can provide theoretical support for the development of grass industry and regional vegetation construction in arid and semi-arid areas.

Identification and analysis of DEGs under drought stress and rehydration
Compared with CK, 2,905 DEGs were identified by padj < 0.05. There were 2,221 DEGs in the whole process of drought stress, among which 1,923 were identified in A vs B (CK vs drought stress 9d), including 1,411 up-regulated genes and 512 down-regulated genes; 489 DEGs were identified in A vs C (CK vs drought stress 12d), including 226 up-regulated genes and 263 down-regulated genes; 924 DEGs were identified in A vs D (CK vs rehydration), including 387 up-regulated genes and 537 down-regulated genes. The results showed that the number of DEGs was the most in A vs B when compared with other 2 treatments. It indicated that expression of genes was abundant and molecular activity was the most active in pathways that response to drought stress during this stage. In A vs C, the number of DEGs of Medicago ruthenica decreased significantly, the plant gradually adapted to the drought stress, and the molecular activity tended to be stable. To sum up, more genes were up-regulated in the leaves of the forage to respond to drought stress, while after rewatering, more genes were down regulated (Fig. 4).
Different treatments also had their own specific DEGs. There were 1563, 270 and 684 DEGs that only expressed in A vs B, A vs C and A vs D respectively (Fig. 5). The results showed that different measures were taken by the plant to respond to distinct external environment.

Go enrichment analysis of DEGs under different treatment
The Blast 2 GO v 2.5 software was used to analyze the GO enrichment of these DEGs in order to well known the biological functions of them. The results showed that 1,319 among 1,923 DEGs were annotated successfully in A vs B, including 987 up-regulated genes and 332 down-regulated genes. In the BP category, the top three enriched terms were "oxidation-reduction process" (126 up-regulated genes, 48 down-regulated genes), "carbohydrate metabolic process" (85 up-regulated genes, 47 down-regulated genes) and "development process" (126 up-regulated genes, 48 down-regulated genes). Most of the DEGs enriched in above biological processes were up-regulated. In addition, DEGs were also significantly enriched in antibiotic, polyol, inositol, glycoside biosynthetic and metabolic process, drug metabolic process and mitochondrial electron transport, cytochrome c to oxygen biological process and they were all up-regulated. Under the "cellular component" group, a lot of DEGs were classified into "respiratory chain complex IV", "respiratory chain", "cytochrome complex", and most of them were up-regulated. Their molecular functions were mainly related to "oxidoreductase activity", "hydrolase activity" and "electron carrier activity". The genes enriched in "oxidoreductase activity" and "electron carrier activity" were mostly up-regulated, while the genes enriched in "hydrolase activity" were mostly down-regulated (Additional file 3: Table S3; Fig. 6). 316 among 489 DEGs were successfully annotated in A vs C, including 132 up-regulated genes and 184 down-regulated genes. At this stage, DEGs were mainly involved in carbohydrate metabolism process (glucan, polysaccharide, sucrose, starch, disaccharide and oligosaccharide metabolic process), and most of the genes were down-regulated. In addition, the DEGs were also significantly enriched in cell wall related biological processes (cell periphery, cell wall organization/modification), and all of them were down-regulated. In the CC category, cell wall and external encapsulating structure were the dominant enriched terms. The molecular functions of these down-regulated DEGs were mainly related to enzyme activities, including pectinesterase activity, hydrolase activity, antioxidant activity, oxidoreductase activity and peroxidase activity etc.. The up-regulated genes of this process were mainly concentrated in inositol metabolic process and a series of catabolic processes, including cellular carbohydrate catabolic process, polyol catabolic process and inositol catabolic process, alcohol catabolic process and organic hydroxy compound catabolic process. And their molecular functions were mainly related to the enzyme activity, including inositol oxygenase activity, glutamate-5-semialdehyde dehydrogenase activity, and glutamate-5-kinase activity and amino acid kinase activity (Additional file 4: Table   S4; Fig. 6). 632 among 924 DEGs were annotated successfully in A vs D, including 263 up-regulated genes and 369 down-regulated genes. The DEGs were mainly involved in single-organism metabolic process, carbohydrate metabolic process, polysaccharide metabolic process, polysaccharide and carbohydrate catabolism, glucan metabolic process, oligosaccharide metabolic process and disaccharide metabolic process. For the "molecular function" group, hydrolase activity, oxidoreductase activity, beta-amylase activity and dioxygenase activity were the most abundant subcategories. The molecular function of the up-regulated genes was mainly related to amylase activity, which were involved in the process of polysaccharide catabolic; the molecular function of the down-regulated genes was related to hydrolase activity, which were involved in carbohydrate metabolic process (Additional file 5: Table S5; Fig. 6).
According to the analysis of the enrichment of KEGG metabolic pathway, unigenes were assigned to 104 pathways and 20 pathways were significantly enriched (p-value < 0.05). Drought  Table S8; Fig. 9).

Quantitative real-time PCR analysis
In order to verify the accuracy of transcriptome sequencing, 6 DEGs were selected for qRT-PCR, including a glucan endo-1,3-beta-glucosidase (Fig. 10A), a ethylene response transcription factor ERF034 (Fig. 10B), a photosynthesis related gene carotene isomerase (Fig. 10C), a lipid metabolism related gene GDSL lipase (Fig. 10D), a DNA damage repair/tolerance protein DRT100 (Fig. 10E), and a fatty acid desaturase (Fig. 10F). The Medicago ruthenica MrActin gene was used as the endogenous reference. The expression patterns detected in qRT-PCR fit well with those in RNA-seq analysis (Fig. 10). Such results demonstrated that DEGs identified based on transcriptome sequencing were reliable.

Disscussion
In this study, we compared the differences of gene expression of Medicago ruthenica seedlings under drought stress (9d, 12d) and rehydration with normal water supply. which showed that the compensation effect of rehydration was remarkable, but it didn't mean 100% repair. Meanwhile, some new different genes were also produced.
According to venn analysis of the DEGs under different treatments, there were 191 DEGs in 9th vs 12th day of drought stress. These genes may be helpful for Medicago ruthenica to cope with different degrees of drought stress, and can be used as candidate genes to improve the stress resistance of other plants. There were 43 DEGs shared by the 3 treatments, which may not only participate in the response of Medicago ruthenica to drought stress, but also joined in the process of restoration after rehydration. These genes can be used as the candidate genes to strengthen the drought adaptability of plants. Moreover, in addition to the common genes, there were also specific DEGs existed in each treatment. 1563, 270 and 684 DEGs were expressed specific on the 9th, 12th day of drought stress and rehydration respectively. The result showed that a great quantity of new DEGs appeared in different treatments, indicating that different strategies may be adopted to adapt the different degrees of drought stress and rehydration. Therefore, we analyzed the adaptive strategies of Medicago ruthenica in different treatment by GO and KEGG.

Enrichment analysis of DEGs
In this study, 2 905 DEGs were screened out and their GO and KEGG enrichment analysis were carried out. The results also proved that Medicago ruthenica adopted distinct strategies to adapt the different degrees of drought stress and rehydration for a further step. After 9 days of drought stress, it was found that the molecular function of DEGs was mainly related to the activities of oxidoreductase, hydrolase and electron transport, and they were mainly enriched in the process of Medicago ruthenica still used the soluble substances produced in the process of polysaccharide catabolism as the sources of energy, carbon and nitrogen. It was convenient for centralized supply of nutrients after rehydration. On the one hand, it was conducive to accelerating the plant recovery process [14]. On the other hand, it played a role of osmotic regulation to strengthen the ability of absorbing water and was helpful to the over recovery of leaf water potential and water transport efficiency.
Under drought stress, plants usually start a complex network regulation mechanism to change the genes' expression level, promote the interaction between biochemical and molecular processes in the body, so as to adapt to the stress environment of water deficiency. In this study, the DEGs and its concentrated metabolic pathways were screened out by transcriptome analysis. It reflected the regulatory process of Medicago ruthenica's response to sustainable drought stress and rehydration on the transcriptional level. It was helpful to preliminarily explained the molecular mechanism of adaptation to drought stress and rehydration. In this period, Medicago ruthenica mainly up regulating the DEGs in "Galactose metabolism", "Arginine and proline metabolism" and "Starch and sucrose metabolism" to play a role in osmoregulation. In this work, α-galactosidase, GOLS, HK and PGM were up-regulated in "Galactose metabolism". Galactose, as a soluble sugar, played an important role in osmoregulation.
The change of gene expression in this pathway could enhance the water holding capacity and improve the drought tolerance of Medicago ruthenica. These genes were in favor of the active adaptation of the forage to arid environment. Simultaneously, 11 P5CS and a OAT were up-regulated in "Arginine and proline metabolism". Proline, as a major osmotic protective substance, was synthesized a lot in arid environment. Its synthesis was connected with the temporary starch degradation, and played a role in plant to cope with drought stress [15]. Proline was the most soluble amino acid with strong hydrophilicity. And it could make for ammonia detoxification and protection of plasma membrane. It was found that there were two pathways for proline synthesis in plants: P5CS and OAT [16,17]. In our study, only one OAT gene was identified in the whole process of drought stress, indicating that the plant mainly synthesized proline through P5CS pathway. So, Medicago ruthenica synthesized a large amount of proline by inducing the up-regulation of P5CS in a water shortage environment. It can not only stabilize the protoplast colloid and metabolism process in tissues, help cells and tissues to retain water, reduce water loss in plants, but also protect the integrity of plasma membrane by detoxification. Besides, β-fructofuranosidase genes, β-glucosidase genes, galacturonidase genes, galU, GPI, HK and PGM in "Starch and sucrose metabolism" were up-regulated in this period. These genes promoted the decomposition of non-structural carbohydrates, produced soluble sugar, protein and released energy when plants suffered from water deficit. The expression of these genes not only played the role of osmotic regulation but also provided energy for plants to resist drought. We also found a up-regulated gene TPS in this pathway. It was a key enzyme gene in trehalose synthesis pathway, and trehalose was related to the tolerance of plants to abiotic stress. Yang [18] found that the length, diameter, number of primary and secondary lateral roots, total root area and root activity of transgenic maize containing TPS1 increased obviously. The expression of this gene might promote the root growth and increase the root/shoot ratio of Medicago ruthenica, which was beneficial to withstand the moderate drought (9d).
During this stage, "Nitrogen metabolism" was inhibited by the down-regulated genes NR and GLT1 and the pathways related to carbohydrate metabolism was promoted in Medicago ruthenica.
It was suggested that the Oxalacetate-Malate Shuttle, which is responsible for the conversion of NADPH produced by light reaction to NADH in cytoplasm, may be a valve on chloroplast membrane. The closed valve inhibited the export of reducing products in chloroplast Malic acid accumulated, CO2 assimilation was blocked and NO2reduction was promoted. Conversely, it was conducive to the assimilation of CO2 and photosynthesis of plants [19]. Therefore, inhibiting Nitrogen metabolism and using more energy to promote carbohydrate metabolism may be one of the mechanisms of Medicago ruthenica to adapt to drought stress and rehydration.
Under the moderate drought stress (9d), the antioxidant defense system was also activated in response to drought stress. In the process of plant growth, drought stress induced plant oxidative stress and produced a lot of ROS. Low concentration of ROS could act as a second messenger to mediate multiple responses in cell signaling pathways, but high concentration of ROS could cause cell damage or even death [20,21]. In order to against ROS toxicity, Medicago ruthenica activated its antioxidant system to maintain or reconstruct redox balance. The excessive ROS in the body was eliminated through various metabolic pathways including Peroxisome, Glutathione metabolism [22]. 3 genes encoding SOD in "Peroxisome" were up-regulated under drought stress. Moreover, Medicago ruthenica also could resist the stress of this stage by slowing down its growth. In our study, photosynthesis of the forage was inhibited under the continuous drought stress, and many related genes were significantly down regulated. This result were similar to those of De [25], Hanf [26] and Behringer [27]. In the forage, a total of 10 LHC genes, including 7 LHCB1 genes and 3 LHCB3 genes, were identified in "Photosynthesis-antenna proteins" under this condition. And all of them were down-regulated. Consequently, the light harvesting ability of Medicago ruthenica was obviously limited and its photosynthesis was suppressed. So, the decrease of photosynthetic capacity of the plant under water deficit was mainly caused by the decrease of light harvesting capacity. The inhibition of "Nitrogen metabolism" also affected the synthesis of chlorophyll to a certain extent. The NADPH-protochlorophyllide oxidoreductase gene in "Porphyrin and chlorophyll metabolism" was down-regulated, which inhibited the production of chlorophyllide. Although the photosynthetic capacity of the plant was inhibited, "Carbon fixation in photosynthetic organisms" was still in progress. There were two MDH, two PGK, a phosphoenolpyruvate carboxykinase and a GAPDH enriched in this pathway and all of them were up-regulated. That is, the photosynthesis of Medicago ruthenica was still in progress so that the plants could grow slowly and the biomass could increase continuously. At the same time, under drought stress, stomata tend to be closed, and the CO2 absorbed from the atmosphere was reduced.
It caused extensive accumulation of NADPH. The synthesis of proline consumed NADPH and ATP and provided NADP + and ADP for photosynthesis of plants. It stimulated the photosynthesis of plant so that Medicago ruthenica could maintain a certain growth rate.
To maintain the above metabolism, the body needs enough energy supply. Under moderate drought stress, energy supply in Medicago ruthenica mainly depends on "TCA cycle" and "Pentose phosphate". Genes in "TCA cycle" were up-regulated, including ACO, IDH, MDH and so on. Meanwhile, genes in "Pentose phosphate", including PGD and PGM, were also up-regulated to generate ATP. It provided energy for the life activities of Medicago ruthenica.
Additionally, INO1, MIP synthase related gene in "Inositol phosphate metabolism", was up-regulated to synthesize myo-inositol. Under the action of TPI, it produced glyceraldehyde-3P and entered the "Glycolysis/gluconeogenesis", which could also provide energy for the plant.
After 12 days of drought stress, most of the genes enriched in "Arginine and proline metabolism", "Cysteine and methionine metabolism", "Ascorbate and aldarate metabolism", "Inositol phosphate metabolism", "Valine, leucine and isoleucine degradation", "Carotenoid biosynthesis" were up-regulated, while those enriched in "Photosynthesis-antenna proteins", "Starch and sucrose metabolism", "Galactose metabolism", "Nitrogen metabolism", "Porphyrin and chlorophyll metabolism", and "TCA cycle", were down-regulated in response to severe drought stress. It proved that serious drought stress suppressed the metabolism of carbohydrate and photosynthesis but promoted the metabolism of amino acid. In conclusion, under severe drought stress, photosynthesis of Medicago ruthenica decreased gradually, and its products decreased accordingly. The forage mainly relied on the synthesis of proline and the decomposition of other stored organic matter to carry out osmotic adjustment, obtain energy and sacrifice part of biomass to maintain physiological metabolism. This was a passive response to drought stress. Apart from that, Medicago ruthenica also synthesized two antioxidants, carotenoids and ascorbic acid, to alleviate oxidative damage and adapt to drought stress. Ultimately, Medicago ruthenica adapted to the drought stress by sacrificing its own biomass.
After rehydration, most of the genes enriched in "Flavonoid biosynthesis", "Circadian rhythm-plant", "Vitamin B6 metabolism", "Carotenoid biosynthesis" and "Diterpene biosynthesis" were up-regulated to participate in the rehydration recovery process of Medicago ruthenica. "Photosynthetic-antenna protein" was restored and the limitation of light capture ability of Medicago ruthenica was relieved. However, "Carbon fixation of photosynthetic tissue" was not completely restored. "Glycerolipid metabolism" was significantly enriched under the condition of rehydration, and the expression of GPAT in it was up-regulated. GPAT was the first enzyme in the pathway of glycerin synthesis, which was helpful for the repair of photosynthetic membrane system and the recovery of photosynthesis. Carbonic anhydrase β-CA5 in N metabolism was significantly down-regulated after rehydration, which maintained the relative stability of photosynthetic rate to a certain extent.
Besides, "Starch and sucrose metabolism" was also not completely recovered. Medicago ruthenica still relied on the decomposition of other nutrients to repair its damage. At this moment, HMGL and hydroxymethylglutaryl CoA synthetase in "Valine, leucine and isoleucine degradation" were up-regulated to generate acetyl coenzyme A. It may be involved in "TCA cycle", "Tryptophan skeleton biosynthesis", "Ketone synthesis and degradation", and could provide C, N and energy for the repair of the herbage. Proline had a unique structure, and could be used as a molecular chaperone to protect the integrity of protein and the activity of enzymes.
When drought stress was relieved, proline in plants could be oxidized and decomposed. This process generates large amounts of FADH2 and NADPH and released a large amount of free energy to generate ATP. It provided the necessary energy for the restoration of plants after stress.
After rehydration, many significantly enriched metabolic pathways under drought stress were recovered, but there were still some metabolic pathways significantly enriched. It is speculated that these metabolic pathways may be beneficial to the recovery process of Medicago ruthenica. adversity [29]. At present, many studies had confirmed that the circadian rhythms of plants were directly related to the response of plants to drought stress [30]. Some critical drought genes and physiological activities induced by drought showed differences between day and night [31], and could regulate the opening and closing of stomata. However, in this study, although the pathway was enriched both in drought stress and rehydration, it was more significantly enriched in the process of rehydration. It was suggested that rehydration had a significant effect on the circadian rhythm of Medicago ruthenica, and the regulation of the pathway may be beneficial to the recovery of its growth. This pathway including up-regulated genes FKF1, PRR5, TOC1 and down-regulated gene LHY. These genes may be involved in the light signal network transmission.
Regulating the circadian rhythm of plants through these differential genes may accelerate the recovery process of Medicago ruthenica. It was speculated that these genes may be the key genes to regulated the circadian rhythm in response to rehydration. The up or down regulation of these genes affected the photosynthetic pathway, flowering, and cell elongation of euphorbiella sinensis.
Moreover, CHS in "Flavonoid biosynthesis" was also involved in this pathway, and its up-regulation may be beneficial to UV-B protection of leaves.  Primer. Finally, PCR products were purified(AMPure XP system) and library quality was assessed (Agilent Bioanalyzer 2100 system).

Sequence assembly and Unigene annotation
Transcriptome assembly was accomplished based on the left.fq and right.fq using Trinity [32] with min_kmer_cov set to 2 by default. All other parameters also set default. Gene function was annotated based on the following databases: NR, NT, Pfam, KOG/COG, Swiss-Prot, KO, GO.

Analysis of DEGs
Using RSEM to estimate gene expression levels [33]. Using the DESeq R package(1. 10.1) to perform the differential expression analysis of different groups. Genes with an adjusted P-value <0.05 found by DESeq were identified as DEGs. GO enrichment analysis of the DEGs was implemented by the GOseq R packages based Wallenius non-central hyper-geometric distribution [34]. KEGG [35] is used for understanding high-level functions and utilities of the biological system from molecular-level information (http://www.genome.jp/kegg/). We used KOBAS [36] software to test the statistical enrichment of DEGs in KEGG pathways.

Real-time quantitative PCR analysis
qRT-PCR analysis was performed to validate the DEGs identified by transcriptome sequencing.
Total RNA were extracted as described before. The qRT-PCR assays were performed using Fast Actin gene was selected as the endogenous reference. All the primers were synthesized by Invitrogen(Beijing, China) ( Table 2). A 20μl reaction mix was set up containing 10 μl 2×Fast Super EvaGreen ® Master Mix, 1μl 10×ROX, 5.5 μl ddH2O, 3 μl total RNA and 0.5 μg of each primer. The relative expression changes of the endogenous reference and tested genes were analyzed by the 2 -△△ CT method [37].