Next Article in Journal
Regulation of von-Willebrand Factor Secretion from Endothelial Cells by the Annexin A2-S100A10 Complex
Next Article in Special Issue
Analysis of the Coding and Non-Coding RNA Transcriptomes in Response to Bell Pepper Chilling
Previous Article in Journal
Recent Advances in Scar Biology
Previous Article in Special Issue
Genome-Wide Identification and Characterization of Warming-Related Genes in Brassica rapa ssp. pekinensis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Whole-Transcriptome Sequence Analysis of Verbena bonariensis in Response to Drought Stress

1
Department of Ornamental Horticulture, Sichuan Agricultural University, 211 Huimin Road, Wenjiang District, Chengdu 611130, Sichuan, China
2
Shanghai Key Laboratory of Plant Functional Genomics and Resources, Shanghai Chenshan Plant Science Research Center, The Chinese Academy of Science, Shanghai Chenshan Botanical Garden, 3888 Huagong Road, Songjiang District, Shanghai 201602, China
*
Authors to whom correspondence should be addressed.
Int. J. Mol. Sci. 2018, 19(6), 1751; https://0-doi-org.brum.beds.ac.uk/10.3390/ijms19061751
Submission received: 20 April 2018 / Revised: 5 June 2018 / Accepted: 8 June 2018 / Published: 13 June 2018
(This article belongs to the Special Issue Plant Genomics)

Abstract

:
Drought is an important abiotic factor that threatens the growth and development of plants. Verbena bonariensis is a widely used landscape plant with a very high ornamental value. We found that Verbena has drought tolerance in production practice, so in order to delve into its mechanism of drought resistance and screen out its drought-resistance genes, we used the RNA-Seq platform to perform a de novo transcriptome assembly to analyze Verbena transcription response to drought stress. By high-throughput sequencing with Illumina Hiseq Xten, a total of 44.59 Gb clean data was obtained from T01 (control group) and T02 (drought experiment group). After assembly, 111,313 unigenes were obtained, and 53,757 of them were annotated by compared databases. In this study, 4829 differentially expressed genes were obtained, of which 4165 were annotated. We performed GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses, and explored a lot of differently expressed genes related to plant energy production, hormone synthesis, cell signal transduction, and metabolism to understand the stress response of Verbena in drought stress. In addition, we also found that a series of TFs related to drought-resistance of Verbena and provide excellent genetic resources for improving the drought tolerance of crops.

1. Introduction

Adverse environmental factors such as low temperature and salt stress, together with drought, prevent plants from realizing their full genetic potential and have been the main problems facing agriculture [1]. Drought stress affects plant growth and development by affecting the plant respiration, growth, photosynthesis, assimilate partitioning, moisture and nutrient relationships, and drought-induced crop yield losses may outweigh losses from all other causes [2]. A series of physiological and biochemical reactions of the plant under drought conditions are altered through gene regulation, such as activation of respiration, repression of cell growth and photosynthesis, and stomatal closure [3]. Accelerating the pace of revealing drought-tolerance mechanisms will greatly help traditional breeding efforts and the application of modern genetic methods in improving the drought tolerance of crops [4].
Many of Verbenaceae’s plants have important medicinal values, such as Aloysia triphylla [5] and Cordia verbenace [6]. Verbena (Verbena bonariensis L.) is an excellent ornamental landscape plant with extensive management. It has a flourishing and long flowering period, such that it plays an extremely important role in landscape layout. In addition to its ornamental value, Verbena has shown good performance under drought stress in production practice, which provides new thinking for our study on plants’ responses to abiotic stress. However, the understanding of its drought-resistance mechanism is still in the early stage. Presently, high-throughput sequencing technology has been widely used to reveal plants’ intrinsic physiological mechanism at the molecular level, such as model plants Nicotiana tabacum L. [7,8,9], Oryza sativa L. [10,11,12], and A. thaliana L. [13], as well as others, like Glycine max (Linn.) Merr. [14], Brassica napus L. [15], and Cucumis sativus L. [16], for its high accuracy and sensitivity of gene discovery. However, in non-model plants, the progress of drought-resistance research has been slow and unevenly developed, and is thus in need of a lot of effort. In this study, we analyzed the transcriptome and differently expressed genes of Verbena by using high-throughput sequencing technology, ultimately analyzing its mechanism of drought resistance and providing potential drought resistance gene information for resistance breeding work.

2. Results

2.1. Phenotypic and Physiological Indicators of Verbena under Drought Stress

In this study, the morphology of Verbena plants was not significantly affected by drought stress (Figure 1a,b). Different from the morphological indicators, physiological indicators of Verbena had undergone significant changes. The chlorophyll content of leaves first decreased rapidly and then increased by a small margin (Figure 2a), the content of proline (Pro) and soluble protein showed a trend of increasing (Figure 2b,c), the content of superoxide dismutase (SOD) reached the highest level on the 10th day, both catalase (CAT) and peroxidase (POD) gradually increased (Figure 2d and Figure 3a,b), malonaldehyde (MDA) content also increased and relative water content (RWC) did not significantly decrease (Figure 3c,d).

2.2. Sequencing and Annotation of Transcription and Unigenes

Based on Sequencing By Synthesis (SBS), six transcriptomes were sequenced by Illumina Hiseq Xten (Illumina, CA, USA). We obtained a total of 44.59 Gb clean data, and in each sample, the Q30 base was not less than 92.87%, the CG (guanine and cytosine basic groups) content was not less than 44.41% (Table SA1). The Pearson’s Correlation Coeffiencient r between T1–T3 and T4–T6 was listed in Figure SA1. Using the de novo assembly program Trinity [17] to assemble short-reads, a total of 258,326 transcripts with an average length of 1139.68 bp were obtained. After continuing to cluster and assemble the transcripts for analysis and a total of 111,313 unigenes with an average length of 697.08 bp and N50 of 1223 bp were obtained, among which 40,340 (36.24%) were over 500 bp in length (Table SA2 and Figure SA2).
To predict and analyze the function of the Verbena unigenes, we used BLAST software to compare the amino acid sequence with NR (NCBI non-redundant), KOG (euKaryotic Orthologous Groups), COG (Clusters of Orthologous Groups), GO (Gene Ontology), KEGG (Kyoto Encyclopedia of Genes and Genomes), Pfam (Protein family) and Swissprot-Annotation (a manually annotated and reviewed protein sequence database) database, setting BLAST parameters E-value ≤ 1 × 105 and HMMER parameters E-value ≤ 1 × 1010 as standard, a total of 53,757 unigenes was obtained, accounting for 48.29% (111,313) of the total. There were 51,352 (95.53%), 28,994 (53.54%), 14,836 (27.60%), 16,938 (31.51%), 20,988 (39.04%), 32,735 (60.89%) and 27,990 (54.53%), unigenes assigned to these databases, respectively (Table SA3). Only 48.29% of unigenes can be matched to known genes, which may be caused by the current lack of studies on Verbena.
A total of 18,104 (35.27%) of unigenes were annotated to Sesamum indicum, which means that Sesamum indicum had the highest level of homology with Verbena, followed by Erythranthe guttata (8.60%) and Erysiphe necator (7.56%). In addition, Verbena and Sesamum indicum are quite similar in morphology because they have spike inflorenscence, which is morphological proof of their high homology (Figure SA3).
There are 28,994 unigenes assigned to KOG database and 14,836 to COG database. In these two databases, the proportion of “Signal transduction mechanisms” related to plant resistance respectively occupied 9.85% and 9.45%. In addition, the amount of unigenes assigned to the classes related to plants’ response to stress, such as “Defense mechanisms”, “Secondary metabolites biosynthesis, transport and catabolism” and “Inorganic ion transport and metabolism”, was 7.40% and 10.50% in the two databases, respectively (Figure SA4).
20,988 unigenes were annotated and classified into 3 categories of GO: cell component (CC), molecular function (MF), and biological processes (BP). Most of the genes were assigned to the biological process (59.16%), followed by the molecular function (23.87%) and cellular component (16.97%) (Figure SA5). Among these, the first three categories of BP were “metabolic process” (2.85%), “cellular process” (2.35%) and “single-organism process” (1.84%).
KEGG is a suite of databases and associated software for understanding and simulating higher-order functional behaviors of cells or the organisms based on their genome information. There were 16,938 (31.51%) unigenes allocated to 129 pathways of KEGG, and the pathways assigned with the most genes are “ribosome” (922), “carbon metabolism” (708) and “biosynthesis of amino acids” (622).

2.3. Analysis of Differentially Expressed Genes (DEGs)

Using Bowtie [18], the clean reads were compared to the unigene library. Genes’ expression levels were estimated by RSEM (RNA-Seq by Expectation Maximization) [19], and according to the results, the expression abundance of a single gene was expressed as the value of FPKM (transcript fragment per million fragments). The volcano map showed that there were 4829 DEGs, where 3841 (79.54%) were up-regulated and 998 (20.46%) were down-regulated (Figure SA6). Finally, 4165 (85.72%) DEGs were annotated, most of which were annotated to NR (4155) and eggNOG (3957), followed by Pfam (3419), Swissprot (2577), KOG (2426), GO (1756), COG (1554) databases and KEGG (1521) (Table SA3). Clustering results of all DEGs are shown in Figure SA7.
Among the DEGs assigned to KOG and COG, the number of genes in these classes related to plant-resistance of abiotic stress such as “inorganic ion transport and metabolism”, “secondary metabolites biosynthesis, transport and catabolism” and “defense mechanisms” were, respectively, 206 (8.50%) and 171 (11.00%), higher than it was for unigenes (Figure 4).
We found that BP, MF and CC accounted for 58.48%, 24.64% and 16.88% respectively by the analysis of 1756 DEGs in GO database. In the entire genetic background and DEGs, the enrichment of genes for each node in the GO database is listed in Figure 5. The top three classes in BP are “DNA integration”, “RNA-dependent DNA replication”, and “photosynthesis, light harvesting”, the most notably of which is the enrichment of “photosynthesis, light harvesting”, and the most noticeable of which in MF are the “acyl-CoA dehydrogenase activity” and “oxidoreductase activity, acting on paired donors, with oxidation of a pair of donors resulting in the reduction of molecular”; the p-values obtained by the KS test were 0.00024 and 0.00016, respectively. “Ribosome”, “fungal-type vacuole membrane” and “nucleus” are the top three in CC, followed up by “photosystem II” and “photosystem I”, with p-values of 8.3 × 105 and 1.4 × 104. The top five significantly enriched nodes of BP, MF and CC are listed in Table SA4.
The DEGs were distributed to 126 lower pathways of KEGG. The first three pathways with the highest enrichment factor were “photosynthesis antenna proteins”, “betalain biosynthesis” and “flavone and flavonol biosynthesis”, with factors of 6.84, 6.36 and 4.71, respectively. The top ten most significantly enriched pathways of DEGs in KEGG are listed in Table SA5 . In the “Photosynthesis antenna proteins” pathway, the class with the highest enrichment factor had a total of 20 DEGs, which were all annotated to the homologous gene annotations of “Light-harvesting Chlorophyll a/b Binding (LHC) Proteins”—the apoproteins of the light-harvesting complex of photosystem II (PSII) [20].

2.4. DEGs of Transcription Factors (TFs) under Drought Stress

In this study, a total of 2146 (44.44%) transcription factor DEGs were identified, of which 1656 (77.17%) were up-regulated and 490 (22.83%) were down-regulated. The TFs mainly focused on bHLH, MYB, ERF, NAC and C2H2 families, with 228 (10.62%), 207 (9.65%), 176 (8.20%), 153 (7.13%) and 102 (4.75%) DEGs (Figure 6), respectively. The number of up-regulated genes of the ten were 171 (75%), 157 (75.85%), 137 (77.84%), 122 (79.74%) and 76 (74.51), respectively.

2.5. Expression Level of DEGs’ Changes and Verification Using qRT-PCR

We selected ten genes involved in different important biological processes to perform qRT-PCR, including genes related to “energy production and conversion”, “transcription factors”, “LHCB2 proteins”, “glutamine synthetase”, “protein kinases”, “lipid transport and metabolism”, “carbohydrate transport and metabolism” and “nitrogen metabolism”. During the experiment, similar to indexes of physiological, the expression level of these DEGs were all up-regulated with the increase of stress time (Figure 7a,b). To further verify the expression level of genes obtained from Illumina Hiseq Xten, we compared the data obtained from the 15th day with sequencing, and the results showed a strong correlation between the two (Figure 8).

3. Discussion

3.1. Morphological and Physiological Index Analysis

Among the ten most abundant nodes in the BP of GO database, the most notable one is “photosynthesis, light harvesting”. Many studies [21,22] have shown that water stress can change the plants’ chlorophyll content, which can indicate the sensitivity of plants to water stress and directly affect the photosynthetic yield. Some studies have shown that drought increases the chlorophyll content of plants’ leaves [23], while some other studies believe it would be gradually decreased [24]. In this study, the chlorophyll content decreased rapidly and then increased by a small margin. Although it is unclear what the mechanism of water stress on chlorophyll content is, the fact is that increasing the chlorophyll content will enhance plants’ endurance to survive in adversity, and it is a kind of ability by which plants can adapt to drought stress, indicating that the Verbena has a certain drought tolerance.
In this experiment, the Pro content, soluble protein content, antioxidant enzyme activity and the MDA content were all increased. Pro and soluble proteins are the most common osmotic pressure regulators in drought-stressed plants. Plants that overproduce Pro and soluble protein might acquire the ability to tolerate environmental stresses such as drought and high salinity [25]. There are a large number of antioxidants that can prevent or repair the damage caused by reactive oxygen species and regulate redox-sensitive signaling pathways, such as POD, SOD, CAT and so on [26]. The degree of membrane damage on the leaves can be demonstrated by the content of MDA, which gradually increased in this study. Under the circumstances of drought force, the more drought-tolerant the plants are, the slower the water content of the leaves decreases. Figure 3d shows that on day 15, the RWC of T02 (stress group) was about 16% lower than that of T01 (control group), which performed better than Hordeum vulgare L. [27], Zea mays L. [28] and Glycine max (Linn.) Merr. [29], indicating that the leaves of Verbena have a certain water retaining capacity.

3.2. The Enrichment and Pathway Analysis of DEGs in GO and KEGG Databases

The chlorophyll content decreased sharply in the early stage of drought stress, and eventually increased to a lower level than T01. The 16 DEGs related to “porphyrin and chlorophyll metabolism” are listed in Table 1, and the enzymes of porphyrin metabolism were almost all up-regulated, while the enzymes of chlorophyll metabolism were almost all down-regulated. Analyzing the genes involved in the regulation of stomatal closure by ABA (abscisic acid), we found that some of the genes’ expression levels of PYL and PP2C were significantly up-regulated (Table 2), and according to the analysis of the “carbon fixation in photosynthetic organisms” pathway, we found that a number of genes related to “C4-Dicarboxylic acid cycle and carbon fixation pathways in prokaryotes” process showed an upward trend. This indicates that under drought stress, the photosynthesis of Verbena may be mainly inhibited by stomatal closure, and at the same time it may provide an adequate carbon source for photosynthesis by enhancing the biological carbon sequestration pathway. The DEGs in “photosystem II” and “photosystem I” of CC also showed a high degree of enrichment, indicating that drought had a great impact on the photosynthesis of Verbena, which is the direction we should focus on. In MF, “acyl-CoA dehydrogenase” and “oxidoreductase” are obviously enriched. Previous study has shown that deficiencies of acyl-CoA dehydrogenases can lead to disorders of fatty acid oxidation, leading to life-threatening metabolic disorders [30]. It is known that drought stress will produce a large amount of reactive oxygen species (ROS) in plants and start the massive production of oxidoreductase such as SOD. Therefore, the increase of “acyl-CoA dehydrogenase” and “oxidoreductase” is of great significance for plants in responding to drought stress.
In the KEGG database, the first three pathways with the highest DEGs enrichment factors are “photosynthesis antenna proteins”, “betalain biosynthesis” and “flavone and flavonol biosynthesis”. After comparison, the homologous genes with almost all the down-regulated DEGs (Table 3) in “photosynthesis antenna proteins” pathway are related to LHC proteins, which may be involved in the drought resistance of plants and play an important role in crop environmental adaptability and yield [31]. Alberte R.S. et al. [32] showed that LHC proteins are a target easily attacked by water stress, and the loss of chlorophyll, the increase of chlorophyll a/b ratio and the decrease of photosynthetic unit under water stress are all caused by the decrease of LHC proteins. Vappaavuovi E. et al. [33] also confirmed that water stress reduced the LHC protein complex. In addition, some studies have showed that LHC proteins may be involved in the partial regulation of ABA signaling and play an active role in guarding cell signaling—these proteins may be induced by ABA and positive regulates ABA to inhibit stomatal opening [34,35]. Therefore, the decrease of LHC proteins will inhibit the production of ROS. However, due to the insufficiency research of LHC proteins, further research on its response to abiotic stress is needed. “Betalain biosynthesis” pathway has three distinctly up-regulated genes (Table 4). Betaine is the trimethyl derivative of the amino acid glycine, an efficient methyl donor that promotes fat metabolism and protein synthesis and the increasing of betaine biosynthesis has been shown to play an important role in osmoregulation of plants, which can help plants to withstand drought stress [36]. The “Flavone and flavonol biosynthesis” pathway has four distinctly up-regulated genes with an enrichment factor of 4.71, and studies have shown that the flavonoid substance may contribute to antioxidant functions in response to drought stress [37]. At present, there is still a lack of research on the changes in the content of secondary metabolites that have economic or medicinal value in Verbena under abiotic stress, but according to the schematic pathway figures of betalain and flavone (Figure 9a,b), we can preliminarily speculate that the content of flavonoids and betaine were up-regulated under drought conditions, while other plants also have similar research conclusions. Xing W. et al. [38] have shown that under drought stress, the endogenous leaf glycine betaine level of A. thaliana L. increased about 18-fold over that in the control plants, and similar results were also found in Pyrus bretschneideri Rehd [39] and Hordeum vulgare L. [40]. From the pathway of “flavone and flavonol biosynthese”, we can see that all DEGs are up-regulated (Table 5) and there are two kinds of useful downstream production—luteolin and quercetin. Studies have shown that drought increased the accumulation of luteolin in Ligustrum lucidum Ait. [37], and quercetin also showed a positive impact on Vigna unguiculata L. Walp. [41], but in Cabernet Sauvignon [42], quercetin showed no obvious change under water stress. Therefore, further experiments are needed to investigate the effects of abiotic stress on secondary metabolite content in Verbena.

3.3. Biological Mechanism of Verbena in Response to Drought Stress

The protective mechanisms of plants in response to stress are regulated by alterations in the expression level of stress-responsive genes. Among the DEGs assigned to KOG and COG, genes related to plants’ stress response accounted for quite a high proportion; 8.50% and 11.00%, respectively. We performed GO and KEGG pathway enrichment analysis and excavated a group of important drought-responsive genes related to multiple biological mechanisms of plant energy production, hormone synthesis, cell signaling, and metabolism. In this study, there are 12, 25 and 24 genes differentially expressed in “glyoxylate and dicarboxylate metabolism”, “citrate cycle (TCA cycle)” and “glycolysis/gluconeogenesis” pathway, respectively, and all of them were up-regulated. This indicates that the energy production of Verbena changed a lot, and the main method of providing ATP was changed from the photosynthetic phosphorylation to the oxidative phosphorylation.
ABA is one of the most important factors in abiotic stress response, and is involved in almost all plant activities, such as photosynthesis, ionic homeostasis, and antioxidant defense [43]. The quantities of gene expression controlling the key enzymes, such as PYL and PP2C in the ABA signaling pathway, have changed a lot. The vicissitudinous DEGs of PYL activate the high expression of PP2C and furthermore promote stomatal closure and decrease transpiration to reduce water loss. Under the “plant hormone signal transduction” pathway, in addition to the enzymes of ABA, the level of genes that control PR1 in the “salicylic acid (SA)”, JAZ and MYC2 in the “jasmonic acid (JA)”, GH3 and SAUR in the “auxin”, and MPK6 together with EBF1_2 in “Ethylene” pathways, were almost all up-regulated (Table 2). This is in line with previous findings that the changes of hormonal synthesis and signal transduction are the conservative mechanisms by which plants response to adverse circumstances. It has been demonstrated that NO (nitric oxide) can act as a signaling molecule to activate ROS-scavenging enzymes in drought stress [44]. Under the pathway of “Nitrogen metabolism”, we observed that the gene expression levels of NR and NirA in the assimilatory nitrate reduction pathway were down-regulated, and the expression of GLT1 and GLUL in the glutamate synthase pathway were up-regulated (Table 6). The glutamate synthase pathway leads to the process of glutamate metabolism and the production of ammonia, and a sufficient amount of oxidized glutamic acid will ensure that there is sufficient carbon skeleton for the tricarboxylic acid cycle to function effectively [45]. Therefore, we can speculate that improving the level of nitrogen metabolism is a very efficient method to help plants resist abiotic stress.
Carbohydrate metabolism is the center of the entire biological metabolism and involves the protein metabolism, lipid metabolism, nucleic acid metabolism and secondary metabolites production. Studies have shown that in rice [46], sugar can be used as a signal to induce the expression of genes associated with abiotic stress. Garg A.K. et al. [47] found that the accumulation of trehalose in rice enabled the transgenic rice to exhibit high salt tolerance, drought tolerance, and low-temperature stress. In this study, there were 28 and 12 DEGs in “Starch and sucrose metabolism” and “Pentose phosphate pathway”, respectively, most of which were up-regulated, which indicates the importance of sugar for Verbena to cope with drought stress. Furthermore, we found a series of pathways for the up-regulation of DEGs expression: protein metabolism such as “Biosynthesis of amino acids”, lipid metabolism such as “Fatty acid metabolism” and “Glycerophospholipid metabolism”, nucleic acid metabolism such as “Pyrimidine metabolism” and “Purine metabolism”, and metabolism of secondary products such as “alpha-Linolenic acid metabolism” and “Flavone and flavonol biosynthesis”. All of these show that carbohydrate metabolism is very important for Verbena to improve its tolerance to adverse environmental conditions.
Moreover, “ubiquitin mediated proteolysis”, an energy-consuming, highly efficient and highly directional protein degradation process, is also noteworthy. It plays an important role in many aspects, such as modulation of the immune and inflammatory responses, the regulation of cell cycle, control of signal transduction pathways, development and differentiation etc. [48,49]. In the present study, there are 16 DEGs in this pathway (Table 7), all of which up-regulate. Half are involved in the regulation of ubiquitin conjugating enzyme, and the others are involved in the regulation of ubiquitin ligase. At present, the understanding of this process is still very limited. The genes participated in this process during drought stress, and the mechanism of the enzymes which were regulated by these genes has yet to be studied.
We found evidence that Verbena responds to drought stress by altering energy synthesis pathways, decreasing transpiration, resetting hormone secretion levels, and increasing cell osmotic pressure and glucose metabolism. In general, Verbena’s defensive response, like most plants under stress conditions, is a process of rebuilding physiologically, biochemically and metabolically, from the growth-oriented to the defensively based one. The genes mentioned in Section 3 are listed in Table 1, Table 2, Table 3, Table 4, Table 5, Table 6 and Table 7 and after counting the TFs to which these genes belong, we constructed the MYB TF phylogenetic tree using the reported abiotic stress genes to prepare for the next work (Figure SA8).

3.4. DEGs of Transcription Factors (TFs) under Drought Stress

There are 228 (10.62%) DEGs that were assigned to bHLH TF. The basic helix-loop-helix proteins are one of the largest transcription factor families and are widely distributed in eukaryotes [50]. Many important drought-tolerant genes have been discovered in the bHLH family. Over-expression of OsbHLH148 in the transgenic rice make the plant more drought-tolerant by regulating the jasmonate signal transduction pathway [51]. Up-regulation of bHLH122 can significantly increase ABA levels in cells and it is a positive regulator of drought, NaCl and osmotic signaling [52]. AtbHLH112 is a nuclear-localized protein induced by salt, drought and ABA, and it can increase proline levels and improve ROS scavenging ability to enhance stress tolerance [53]. At the same time, there are also a number of other TFs that play an important role in plants’ resistance to abiotic stresses, such as MYB (9.65%), ERF (8.20%) and NAC (7.13%). In this study, most of the genes that enhance plants’ stress resistance were achieved by overexpression, so we will pay more attention to the up-regulated genes and dig out the role of specific genes in the processes of stress resistance by transgenic and gene-silencing technology.

4. Materials and Methods

4.1. Plant Materials and Drought Treatments

Verbena bonariensis L. with the age of one month were used in this study. The experimental materials were incubated in greenhouse (25 °C, 16 h photoperiod, 50% RH), soil texture for the medium loam, the maximum soil moisture content of 80%. A total of 90 strains of Verbena were randomly divided into two groups T01 (control group) and T02 (drought experiment group), 1 seeding in each pot and 15 in one duplicated group, each of the T01 and T02 contained 3 duplicated groups named T1–T3 and T4–T6, respectively. Determination of soil moisture content was measured by the oven drying method [54]. The soil water content of T01 was kept at 80% of saturated soil moisture all the time, while watering of the T02 was stopped until the soil water content had reduced to 25% of the soil saturated water content. After that, the soil moisture content was measured every 2 days to replenish the amount of deficiency and maintained this drought status for 15 days.
At the 5th, 10th and 15th day, the mature leaves (3rd to 8th functional leaf) were selected randomly from the plants of T01 and T02 for the determination of physiological indexes. After the drought had been maintained for 15 days, samples were collected, rapidly frozen with liquid nitrogen, stored at −80 °C, and finally sent to Biomarker Technologies Co., Ltd. (Beijing, China) for whole-transcriptome sequencing.

4.2. Determination of Morphological and Physiological Characters

Vernier calipers were used to measure root length (whichever is longer) and stem length. Chlorophyll content was determined by NanoPhotometer® spectrophotometer (Implen, CA, USA) [55]; Pro content was measured using acidic-ninhydrin method [56]; estimation of soluble protein by Bradford method [57]; nitroblue tetrazolium blue (NBT) reduction method was used to determine the activity of SOD [58]; CAT activity was determined using UV absorption method [59]; POD activity was determined by guaiacol method [60]; determination of MDA content by thiobarbituric acid test (TBA) [61]. Relative water content (RWC) of leaves was measured according to the method of Tambussi E.A. et al. [62] using the formula: RWC = (FW − DW)/(turgid weight − DW) × 100%, where FW is the leaf fresh weight, DW is the leaf dry weight at 85 °C for 3 d, and SW is the turgid weight of leaves after soaking in water for 4 h at room temperature (approximately 25 °C). Physiological measurements were set up for three replicates to reduce the error.

4.3. Extraction of RNA, Library Preparation for Transcriptome Sequencing

To ensure the qualified samples were obtained for transcriptome sequencing, total RNA was extracted with Trizol kit (Invitrogen, Carlsbad, CA, USA) and its purity concentration and integrity detected by the Nanodrop, Qubit 2.0, Agilent 2100 method. Then, a total amount of 3 μg qualified RNA per sample was used as input material for the RNA sample preparations. According to manufacturer’s recommendations, sequencing libraries were generated using NEBNext®Ultra™ RNA Library Prep Kit for Illumina® (New England Biolabs, Ipswich, MA, USA). The messenger RNAs (mRNAs) were separated from the total RNA by Oligo (dT) and were cleaved into short fragments at random. The first strand cDNA was synthesized by random hexamer primer, then the buffer, dNTPs, DNA polymerase I and RNase H were used to synthesize the second strand cDNAs. Lastly, the cDNAs were purified with AMPure XP beads and after end-repair and single nucleotide A (adenine) addition, the qualified cDNA libraries were constructed by PCR enrichment. After the cDNA libraries were constructed, Qubit 2.0 was used for preliminary quantification, and then the Agilent 2100 was used to detect the insert size of the libraries. After that, the Q-PCR method was used to accurately quantify the effective concentration of the libraries (effective library concentration > 2 nM) to ensure library quality. After passing the screening, high-throughput sequencing was performed with Illumina Hiseq Xten. The raw sequencing data have been submitted to the National Center for Biotechnology Information Search database (NCBI) Sequence Read Archive database with accession number SRP132610.

4.4. Transcriptome Assembly and Gene Functional Annotation

A total of 44.59 Gb of clean data was removed from the original data, including low-quality data. Then Trinity software was used to break reads into shorter K-mers, extend them to obtain Contigs, gather a collection of Contig clusters, and finally obtain the transcript sequence by using de Bruijin graphs algorithm method and reads.
Using BLAST software, the sequencing of the Unigene sequence was compared with NR, KOG, COG, GO, KEGG, Pfam and Swissprot-Annotation databases. After the prediction of the amino acid sequence of the unigene, HMMER software [63] was used to compare it with the Pfam database [64] to obtain functional annotation information of all unigenes.

4.5. Differential Expression Analysis

Differential expression analysis of the sample groups was performed using DESeq [65] to obtain a set of DEGs between the drought group and the control group, and the False Discovery Rate (FDR) was used to correct the p-value of the multiple-hypothesis test. In this study, FDR < 0.01 and FC (Fold Change) ≥ 4 were set as the threshold for significantly differential expression, and could be used to prepare for the follow-up study on the drought stress of Verbena.

4.6. Quantitative Real-Time PCR Analysis

In order to clarify the stress response variation tendency of Verbena at the transcriptional level on different stages of drought stress, we selected 10 DEGs to participate in different important biological processes for qRT-PCR with SsoFast™ EvaGreen® Supermix every 5 days. A 20 μL fluorescent quantitative reaction system contains 10 μL of stain, 2 μL of cDNA template and 300 nM of primers. The PCR settings are as follows:
TemperatureTimeCycle
95 °C30 s
95 °C15 s40 cycles
60 °C30 s
Relative expression levels were calculated by the 2−ΔΔCt method, and a β-actin gene of Verbena bonariensis (Forward primer: GAAAGATGGCTGGAAGAGGG, Reverse primer: GCTATGAACTCCCTGATGGTC) was used as the reference for quantitative expression analysis. The expression pattern of DEGs was analyzed by melting furnace curve. The qRT-PCR primers are shown in Table SA6.

5. Conclusions

The development of plant genomics plays an important role in the effective use of modern molecular biology methods for the genetic improvement of species. After sequencing the transcriptome of drought-stressed Verbena materials by high-throughput sequencing and verifying the results by qRT-PCR, we finally identified DEGs and TFs related to stress resistance and analyzed their biomodulation mechanism, and this proved Verbena’s commendable drought tolerance at physiological and molecular levels. These data will provide excellent genetic resources for improving the drought tolerance of crops and lay a good foundation for the follow-up study of Verbena.

Supplementary Materials

Author Contributions

B.W., X.-Q.L., and Q.-L.L. conceived and designed the experiments. B.W., X.-Q.L., M.-S.X. and Q.-L.L. performed the experiments. B.W., F.-L.L., Q.Z., L.H., Q.-L.L., L.Z., Y.J. and F.Z. analyzed the data. Q.-L.L. contributed reagents/ materials/ analysis tools. B.W. wrote the paper.

Funding

This research was funded by National Natural Science Foundation of China, grant number (31770742) and Innovative Training Program of Sichuan Agricultural Uniersity, grant number (201710626086).

Acknowledgments

This study was supported by Sichuan Agricultural University Ornamental Horticulture lab. We would like to acknowledge the contribution of Qing-Lin Liu for the provision of experimental materials and instruments.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zhu, J.K. Salt and drought stress signal transduction in plants. Annu. Rev. Plant Biol. 2002, 53, 247–273. [Google Scholar] [CrossRef] [PubMed]
  2. Farooq, M.; Wahid, A.; Kobayashi, N.; Fujita, D.; Basra, S.M.A. Plant drought stress: Effects, mechanisms and management. Agron. Sustain. Dev. 2009, 29, 185–212. [Google Scholar] [CrossRef]
  3. Shinozaki, K.; Yamaguchi-Shinozaki, K. Gene networks involved in drought stress response and tolerance. J. Exp. Bot. 2007, 58, 221–227. [Google Scholar] [CrossRef] [PubMed]
  4. Xiong, L.; Wang, R.; Mao, G.; Koczan, J.M. Identification of drought tolerance determinants by genetic analysis of root response to drought stress and abscisic acid. Plant Physiol. 2006, 142, 1065–1074. [Google Scholar] [CrossRef] [PubMed]
  5. Mengistu, F.; Solomon, T.; Tsion, T.; Diriba, G.; Helen, G. In vitro protocol optimization for micropropagation of elite lemmon verbena (aloysia triphylla). Afr. J. Plant Sci. 2017, 11, 369–376. [Google Scholar] [CrossRef]
  6. Sertié, J.A.; Basile, A.C.; Panizza, S.; Matida, A.K.; Zelnik, R. Pharmacological assay of cordia verbenacea; part 1. anti-inflammatory activity and toxicity of the crude extract of the leaves. Planta Med. 1988, 54, 7–10. [Google Scholar] [CrossRef] [PubMed]
  7. Ni, L.; Ren, X.; Xiang, Z.; Wan, W.; Yang, D. Sequencing and characterization of leaf transcriptomes of six diploid Nicotiana species. J. Biol. Res. 2016, 23, 1–12. [Google Scholar] [CrossRef]
  8. Xing, X.; Li, X.; Zhang, M.; Wang, Y.; Liu, B.; Xi, Q.; Zhao, K.; Wu, Y.; Yang, T. Transcriptome analysis of resistant and susceptible tobacco (Nicotiana tabacum) in response to root-knot nematode Meloidogyne incognita infection. Biochem. Biophys. Res. Commun. 2016, 482, 1114–1121. [Google Scholar] [CrossRef] [PubMed]
  9. Bokvaj, P.; Hafidh, S.; Honys, D. Transcriptome profiling of male gametophyte development in Nicotiana tabacum. Genom. Data 2014, 3, 106–111. [Google Scholar] [CrossRef] [PubMed]
  10. Kim, S.; Park, J.; Lee, J.; Shin, D.; Park, D.S.; Lim, J.S.; Choi, I.Y.; Seo, Y.S. Understanding pathogenic Burkholderia glumae metabolic and signaling pathways within rice tissues through in vivo transcriptome analyses. Gene 2014, 547, 77–85. [Google Scholar] [CrossRef] [PubMed]
  11. Mohanty, B.; Kitazumi, A.; Cheung, C.Y.M.; Lakshmanan, M.; de Los Reyes, B.G.; Jang, I.C.; Lee, D.Y. Identification of candidate network hubs involved in metabolic adjustments of rice under drought stress by integrating transcriptome data and genome-scale metabolic network. Plant Sci. 2016, 242, 224–239. [Google Scholar] [CrossRef] [PubMed]
  12. Li, X.; He, Y.; Yang, J.; Jia, Y.H.; Zeng, H.L. Gene mapping and transcriptome profiling of a practical photo-thermo-sensitive rice male sterile line with seedling-specific green-revertible albino leaf. Plant Sci. 2018, 266, 37–45. [Google Scholar] [CrossRef] [PubMed]
  13. Boudichevskaia, A.; Heckwolf, M.; Althaus, L.; Kaldenhoff, R. Transcriptome analysis of the aquaporin AtPIP1;2 deficient line in Arabidopsis thaliana. Genom. Data 2015, 4, 162–164. [Google Scholar] [CrossRef] [PubMed]
  14. Prince, S.J.; Joshi, T.; Mutava, R.N.; Syed, N.; Joao Vitor, M.S.; Patil, G.; Song, L.; Wang, J.; Lin, L.; Chen, W.; et al. Comparative analysis of the drought-responsive transcriptome in soybean lines contrasting for canopy wilting. Plant Sci. 2015, 240, 65–78. [Google Scholar] [CrossRef] [PubMed]
  15. Kubala, S.; Garnczarska, M.; Wojtyla, Ł.; Clippe, A.; Kosmala, A.; Żmieńko, A.; Lutts, S.; Quinet, M. Deciphering priming-induced improvement of rapeseed (Brassica napus L.) germination through an integrated transcriptomic and proteomic approach. Plant Sci. 2015, 231, 94–113. [Google Scholar] [CrossRef] [PubMed]
  16. Pawełkowicz, M.; Zieliński, K.; Zielińska, D.; Pląder, W.; Yagi, K.; Wojcieszek, M.; Siedlecka, E.; Bartoszewski, G.; Skarzyńska, A.; Przybecki, Z. Next generation sequencing and omics in cucumber (Cucumis sativus L.) breeding directed research. Plant Sci. 2015, 242, 77–88. [Google Scholar] [CrossRef]
  17. Grabherr, M.G.; Grabherr, M.G.; Haas, B.J.; Yassour, M.; Levin, J.Z.; Thompson, D.A.; Amit, I.; Adiconis, X.; Fan, L.; Raychowdhury, R.; et al. Full length transcriptome assembly from RNA Seq data without a reference genome. Nat. Biotechnol. 2011, 29, 644–652. [Google Scholar] [CrossRef] [PubMed]
  18. Langmead, B.; Trapnell, C.; Pop, M.; Salzberg, S.L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10, R25. [Google Scholar] [CrossRef] [PubMed]
  19. Li, B.; Colin, N.D. RSEM: Accurate transcript quantification from RNA Seq data with or without a reference genome. BMC Bioinf. 2011, 12, 323. [Google Scholar] [CrossRef] [PubMed]
  20. Liu, R.; Xu, Y.H.; Jiang, S.C.; Lu, K.; Lu, Y.F.; Feng, X.J.; Wu, Z.; Liang, S.; Yu, Y.T.; Wang, X.F.; et al. Light-harvesting chlorophyll a/b-binding proteins, positively involved in abscisic acid signalling, require a transcription repressor, WRKY40, to balance their function. J. Exp. Bot. 2013, 64, 5443–5456. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  21. Brix, H. The Effect of Water Stress on the Rates of Photosynthesis and Respiration in Tomato Plants and Loblolly Pine Seedlings. Physiol. Plant. 1962, 15, 10–20. [Google Scholar] [CrossRef]
  22. Pelleschi, S.; Rocher, J.P.; Prioul, J.L. Effect of water restriction on carbohydrate metabolism and photosynthesis in mature maize leaves. Plant Cell Environ. 1997, 20, 493–503. [Google Scholar] [CrossRef] [Green Version]
  23. Jackson, R.D.; Idso, S.B.; Reginato, R.J.; Pinter, P.J. Canopy temperature as a crop water stress indicator. Water Resour. Res. 1981, 17, 1133–1138. [Google Scholar] [CrossRef]
  24. Huang, B.; Fry, J.; Wang, B. Water relations and canopy characteristics of tall fescue cultivars during and after drought stress. HortScience 1998, 33, 245–256. [Google Scholar]
  25. Wang, B.M.; Chen, J.J.; Chen, L.S.; Wang, X.N.; Wang, R.; Ma, L.; Peng, S.F.; Luo, J.; Chen, Y.Z. Combined drought and heat stress in camellia oleifera, cultivars: Leaf characteristics, soluble sugar and protein contents, and rubisco gene expression. Trees 2015, 29, 1483–1492. [Google Scholar] [CrossRef]
  26. Weydert, C.J.; Cullen, J.J. Measurement of superoxide dismutase, catalase and glutathione peroxidase in cultured cells and tissue. Nat. Protoc. 2010, 5, 51–66. [Google Scholar] [CrossRef] [PubMed]
  27. Teulat, B.; Monneveux, P.; Wery, J.; Borries, C.; Souyris, I.; Charrier, A.; This, D. Relationships between relative water content and growth parameters under water stress in barley: A QTL study. New Phytol. 1997, 137, 99–107. [Google Scholar] [CrossRef]
  28. Zygielbaum, A.I.; Gitelson, A.A.; Arkebauer, T.J.; Rundquist, D.C. Non-destructive detection of water stress and estimation of relative water content in maize. Geophys. Res. Lett. 2009, 36, 91–100. [Google Scholar] [CrossRef]
  29. Singh, B.B.; Gupta, D.P. Proline accumulation and relative water content in soya bean (glycine max) varieties under water stress. Ann. Bot. 1983, 52, 109–110. [Google Scholar] [CrossRef]
  30. Thorpe, C.; Kim, J.J. Structure and mechanism of action of the acyl-CoA dehydrogenases. FASEB J. 1995, 9, 718–725. [Google Scholar] [CrossRef] [PubMed]
  31. Ganeteg, U.; Külheim, C.; Andersson, J.; Jansson, S. Is each light-harvesting complex protein important for plant fitness? Plant Physiol. 2004, 134, 502–509. [Google Scholar] [CrossRef] [PubMed]
  32. Alberte, R.S.; Thornber, J.P. Water stress effects on the content and organization of chlorophyll in Mesphyll and bundle sheath chloroplasts of maize. Plant Physiol. 1997, 59, 351–353. [Google Scholar] [CrossRef]
  33. Vapaavuori, E.; Nurmi, A. Chlorophyll-protein complexes in Salix sp. “aquatica gigantean” under strong and weak light. II. Effect of water stress on the chlorophyll-protein complexes and chloroplast ultrastructure. Plant Cell Physiol. 1982, 23, 791–801. [Google Scholar] [CrossRef]
  34. Zhang, X.; Wang, H.; Takemiya, A.; Song, C.P.; Kinoshita, T.; Shimazaki, K. Inhibition of blue light-dependent H+ pumping by abscisic acid through hydrogen peroxide-induced dephosphorylation of the plasma membrane H+-ATPase in guard cell protoplasts. Plant Physiol. 2004, 136, 4150–4158. [Google Scholar] [CrossRef] [PubMed]
  35. Xu, Y.H.; Liu, R.; Yan, L.; Liu, Z.Q.; Jiang, S.C.; Shen, Y.Y.; Wang, X.F.; Zhang, D.P. Light-harvesting chlorophyll a/b-binding proteins are required for stomatal response to abscisic acid inArabidopsis. J. Exp. Bot. 2012, 63, 1095–1106. [Google Scholar] [CrossRef] [PubMed]
  36. Ratriyanto, A.; Mosenthin, R.; Bauer, E.; Eklund, M. Metabolic, osmoregulatory and nutritional functions of betaine in monogastric animals. Asian-Aust. J. Anim. Sci. 2009, 22, 1461–1476. [Google Scholar] [CrossRef]
  37. Tattini, M.; Galardi, C.; Pinelli, P.; Massai, R.; Remorini, D.; Agati, G. Differential accumulation of flavonoids and hydroxycinnamates in leaves of Ligustrum vulgare under excess light and drought stress. New Phytol. 2004, 163, 547–561. [Google Scholar] [CrossRef] [Green Version]
  38. Xing, W.; Rajashekar, C.B. Glycine betaine involvement in freezing tolerance and water stress in arabidopsis thaliana. Environ. Exp. Bot. 2001, 46, 21–28. [Google Scholar] [CrossRef]
  39. Gao, X.P.; Yan, J.Y.; Liu, E.K.; Shen, Y.Y.; Lu, Y.F.; Zhang, D.P. Water stress induces in pear leaves the rise of betaine level that is associated with drought tolerance in pear. J. Hortic. Sci. Biotechnol. 2004, 79, 114–118. [Google Scholar] [CrossRef]
  40. Hitz, W.D.; Ladyman, J.A.R.; Hanson, A.D. Betaine Synthesis and Accumulation in Barley during Field Water-Stress. Crop Sci. 1982, 22, 47–54. [Google Scholar] [CrossRef]
  41. Goufo, P.; Moutinhopereira, J.M.; Jorge, T.F.; Correia, C.M.; Oliveira, M.R.; Eas, R.; António, C.; Trindade, H. Cowpea (Vigna Unguiculatal. walp.) metabolomics: Osmoprotection as a physiological strategy for drought stress resistance and improved yield. Front. Plant Sci. 2017, 8, 586. [Google Scholar] [CrossRef] [PubMed]
  42. Quiroga, A.M. Water stress and abscisic acid exogenous supply produce differential enhancements in the concentration of selected phenolic compounds in cabernet sauvignon. J. Berry Res. 2015, 2, 33–44. [Google Scholar] [CrossRef]
  43. Cutler, S.R.; Rodriguez, P.L.; Finkelstein, R.R.; Abrams, S.R. Abscisic acid: Emergence of a core signaling network. Annu. Rev. Plant Biol. 2010, 61, 651–679. [Google Scholar] [CrossRef] [PubMed]
  44. Qiao, W.; Li, C.; Fan, L.M. Cross-talk between nitric oxide and hydrogen peroxide in plant responses to abiotic stresses. Environ. Exp. Bot. 2014, 100, 84–93. [Google Scholar] [CrossRef]
  45. Robinson, S.A.; Slade, A.P.; Fox, G.G.; Phillips, R.; Ratcliffe, R.G.; Stewart, G.R. The role of glutamate dehydrogenase in plant nitrogen metabolism. Plant Physiol. 1991, 95, 509–516. [Google Scholar] [CrossRef] [PubMed]
  46. Joo, J.; Lee, Y.H.; Kim, Y.; Nahm, B.H.; Song, S.I. Abiotic stress responsive rice ASR1 and ASR3 exhibit different tissue-dependent sugar and hormone-sensitivities. Mol. Cells 2013, 35, 421–435. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  47. Garg, A.K.; Kim, J.K.; Owens, T.J.; Ranwala, A.P.; Choi, Y.D.; Kochian, L.V.; Wu, R.J. Trehalose accumulation in rice plants confers high tolerance levels to different abiotic stresses. Proc. Natl. Acad. Sci. USA 2002, 99, 15898–15903. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Ciechanover, A.; Orian, A.; Schwartz, A.L. Ubiquitin-mediated proteolysis: Biological regulation via destruction. BioEssays 2000, 22, 442–451. [Google Scholar] [CrossRef]
  49. Olson, B.L.; Hock, M.B.; Ekholm, R.S.; Wohlschlegel, J.A.; Dev, K.K.; Kralli, A.; Reed, S.I. Scfcdc4 acts antagonistically to the pgc-1alpha transcriptional coactivator by targeting it for ubiquitin-mediated proteolysis. Genes Dev. 2008, 22, 252–264. [Google Scholar] [CrossRef] [PubMed]
  50. Ji, X.; Nie, X.; Liu, Y.; Zheng, L.; Zhao, H.; Zhang, B.; Huo, L.; Wang, Y. A bHLH gene from Tamarix hispida improves abiotic stress tolerance by enhancing osmotic potential and decreasing reactive oxygen species accumulation. Tree Physiol. 2016, 36, 193–207. [Google Scholar] [CrossRef] [PubMed]
  51. Seo, J.S.; Joo, J.; Kim, M.J.; Kim, Y.K.; Nahm, B.H.; Song, S.I.; Cheong, J.J.; Lee, J.S.; Kim, J.K.; Choi, Y.D. OsbHLH148, a basic helix-loop-helix protein, interacts with OsJAZ proteins in a jasmonate signaling pathway leading to drought tolerance in rice. Plant J. Cell Mol. Biol. 2011, 65, 907–921. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  52. Liu, W.; Tai, H.; Li, S.; Gao, W.; Zhao, M.; Xie, C.; Li, W.X. bHLH122, is important for drought and osmotic stress resistance in Arabidopsis, and in the repression of ABA catabolism. New Phytol. 2014, 201, 1192–1204. [Google Scholar] [CrossRef] [PubMed]
  53. Liu, Y.; Ji, X.; Nie, X.; Qu, M.; Zheng, L.; Tan, Z.; Zhao, H.; Huo, L.; Liu, S.; Zhang, B.; et al. Arabidopsis AtbHLH112 regulates the expression of genes involved in abiotic stress tolerance by binding to their E-box and GCG-box motifs. New Phytol. 2015, 207, 692–709. [Google Scholar] [CrossRef] [PubMed]
  54. O’Kelly, B.C. Accurate Determination of Moisture Content of Organic Soils Using the Oven Drying Method. Dry. Technol. 2004, 22, 1767–1776. [Google Scholar] [CrossRef] [Green Version]
  55. Chen, Y.E.; Liu, W.J.; Su, Y.Q.; Cui, J.M.; Zhang, Z.W.; Yuan, M.; Zhang, H.Y.; Yuan, S. Different response of photosystem II to short and long-term drought stress in Arabidopsis thaliana. Physiol. Plant. 2016, 158, 225–235. [Google Scholar] [CrossRef] [PubMed]
  56. Zhi, M.; Li, X. Improvement on the method for measuring proline content. Plant Physiol. Commun. 2005, 41, 355–357. [Google Scholar]
  57. Bradford, M.M. A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding. Anal. Biochem. 1976, 72, 248–254. [Google Scholar] [CrossRef]
  58. Durak, I.; Yurtarslanl, Z.; Canbolat, O.; Akyol, O. A methodological approach to superoxide dismutase (SOD) activity assay based on inhibition of nitroblue tetrazolium (NBT) reduction. Clin. Chim. Acta 1993, 214, 103–104. [Google Scholar] [CrossRef]
  59. Zhang, L.; Zhang, L.; Xi, D.; Luo, L.; Meng, F.; Li, Y.; Wu, C.A.; Guo, X. Cotton GhMPK2 is involved in multiple signaling pathways and mediates defense responses to pathogen infection and oxidative stress. FEBS J. 2011, 278, 1367–1378. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  60. Ranieri, A.; Petacco, F.; Castagna, A.; Soldatini, G.F. edox state and peroxidase system in sunflower plants exposed to ozone. Plant Sci. 2009, 159, 159–167. [Google Scholar] [CrossRef]
  61. Ohya, T. Reactivity of alkanals towards malondialdehyde (MDA) and the effect of alkanals on MDA determination with a thiobarbituric acid test. Biol. Pharm. Bull. 1993, 16, 1078–1082. [Google Scholar] [CrossRef] [PubMed]
  62. Tambussi, E.A.; Nogués, S.; Araus, J.L. Ear of durum wheat under water stress: Water relations and photosynthetic metabolism. Planta 2005, 221, 446–458. [Google Scholar] [CrossRef] [PubMed]
  63. Eddy, S.R. Profile hidden Markov models. Bioinformatics 1998, 14, 755–763. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  64. Finn, R.D.; Bateman, A.; Clements, J.; Coggill, P.; Eberhardt, R.Y.; Eddy, S.R.; Heger, A.; Hetherington, K.; Holm, L.; Mistry, J. Pfam: The protein families database. Nucl. Acids Res. 2014, 42, D222. [Google Scholar] [CrossRef] [PubMed]
  65. Anders, S.; Huber, W. Differential expression analysis for sequence count data. Genome Biol. 2010, 11, R106. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Morphological indexes of Verbena under drought stress. T01 is the control group and T02 is the drought experiment group: (a) The determination of stem length; (b) The determination of root length.
Figure 1. Morphological indexes of Verbena under drought stress. T01 is the control group and T02 is the drought experiment group: (a) The determination of stem length; (b) The determination of root length.
Ijms 19 01751 g001
Figure 2. Physiological indexes of Verbena under drought stress. T01 is the control group and T02 is the drought experiment group: (a) The content of chlorophyll; (b) The content of Pro; (c) The content of soluble protein; (d) The content of SOD.
Figure 2. Physiological indexes of Verbena under drought stress. T01 is the control group and T02 is the drought experiment group: (a) The content of chlorophyll; (b) The content of Pro; (c) The content of soluble protein; (d) The content of SOD.
Ijms 19 01751 g002
Figure 3. Physiological indexes of Verbena under drought stress. T01 is the control group and T02 is the drought experiment group: (a) The activity of CAT; (b) The activity of POD; (c) The activity of MDA; (d) The relative water content of leaves.
Figure 3. Physiological indexes of Verbena under drought stress. T01 is the control group and T02 is the drought experiment group: (a) The activity of CAT; (b) The activity of POD; (c) The activity of MDA; (d) The relative water content of leaves.
Ijms 19 01751 g003aIjms 19 01751 g003b
Figure 4. The bar chart of unigenes’ functional classification annotated in COG and KOG databases. The abscissa is the function classifications of COG and KOG databases and the ordinate is the number of DEGs annotated in it.
Figure 4. The bar chart of unigenes’ functional classification annotated in COG and KOG databases. The abscissa is the function classifications of COG and KOG databases and the ordinate is the number of DEGs annotated in it.
Ijms 19 01751 g004
Figure 5. The bar chart of DEGs annotated in the GO classification. The ordinate at the left represents the percentage of the number of genes, the right ordinate represents the number of genes. The above of two ordinates is the number of DEGs, the following is the number of all genes. The abscissa is the classification of GO. The dark bar represents the number and proportion of DEGs that are enriched in GO function, and the light bar represents the number and proportion of genes that are enriched for each GO function.
Figure 5. The bar chart of DEGs annotated in the GO classification. The ordinate at the left represents the percentage of the number of genes, the right ordinate represents the number of genes. The above of two ordinates is the number of DEGs, the following is the number of all genes. The abscissa is the classification of GO. The dark bar represents the number and proportion of DEGs that are enriched in GO function, and the light bar represents the number and proportion of genes that are enriched for each GO function.
Ijms 19 01751 g005
Figure 6. The sector diagram of TFs’ classification and number in DEGs in response to drought stress.
Figure 6. The sector diagram of TFs’ classification and number in DEGs in response to drought stress.
Ijms 19 01751 g006
Figure 7. The line chart of the expression level of the 10 DEGs varied with the degree of drought during the experiment. (a) Five DEGs varied with the degree of drought during the experiment; (b) another five DEGs varied with the degree of drought during the experiment.
Figure 7. The line chart of the expression level of the 10 DEGs varied with the degree of drought during the experiment. (a) Five DEGs varied with the degree of drought during the experiment; (b) another five DEGs varied with the degree of drought during the experiment.
Ijms 19 01751 g007
Figure 8. The bar chart of results of qRT-PCR in 15th day. The relative expression level of ten DEGs identified in the comparison between RNA-Seq and qRT-PCR. The genes relative expression level were determined by 2ΔΔCt method.
Figure 8. The bar chart of results of qRT-PCR in 15th day. The relative expression level of ten DEGs identified in the comparison between RNA-Seq and qRT-PCR. The genes relative expression level were determined by 2ΔΔCt method.
Ijms 19 01751 g008
Figure 9. Two schematic pathway figures of “Betalain biosynthesis” and “Flavone and flavonol biosynthesis”: (a) Effect of drought stress on the expression of genes associated with betaine. (b) Luteolin and quercetin metabolism.
Figure 9. Two schematic pathway figures of “Betalain biosynthesis” and “Flavone and flavonol biosynthesis”: (a) Effect of drought stress on the expression of genes associated with betaine. (b) Luteolin and quercetin metabolism.
Ijms 19 01751 g009
Table 1. DEGs (analysis of differentially expressed genes) in “Porphyrin and chlorophyll metabolism” pathway of KEGG (Kyoto Encyclopedia of Genes and Genomes).
Table 1. DEGs (analysis of differentially expressed genes) in “Porphyrin and chlorophyll metabolism” pathway of KEGG (Kyoto Encyclopedia of Genes and Genomes).
TermGene IDlog2FCGene DescriptionFDR
URODc76376.graph_c02.33uroporphyrinogen decarboxylase chloroplast precursor1.44 × 10−5
COX15c59080.graph_c02.40uroporphyrinogen decarboxylase chloroplast precursor2.71 × 10−8
FECHc69481.graph_c02.58protoporphyrin/coproporphyrin ferrochelatase7.56 × 10−5
c86128.graph_c23.02chloroplastic isoform X21.01 × 10−7
EARSc69469.graph_c02.43glutamyl-tRNA reductase1.14 × 10−6
c72758.graph_c02.31Porphyrin and chlorophyll metabolism2.72 × 10−3
hemAc85183.graph_c04.04glutamyl-tRNA reductase 1, chloroplastic-like1.98 × 10−11
c77400.graph_c0−2.69hypothetical protein9.98 × 10−38
c77400.graph_c1−2.65glutamyl-tRNA reductase 1, chloroplastic1.84 × 10−31
c77400.graph_c2−2.68glutamyl-tRNA reductase 1, chloroplastic1.67 × 10−23
chlHc88820.graph_c1−2.50magnesium chelatase subunit H1.91 × 10−60
chlEc77176.graph_c0−2.24magnesium-protoporphyrin IX monomethyl ester (oxidative) cyclase7.54 × 10−37
porc85861.graph_c0−3.48protochlorophyllide reductase1.36 × 10−11
chlPc80298.graph_c1−2.94geranylgeranyl diphosphate/geranylgeranyl-bacteriochlorophyllide a reductase1.13 × 10−47
Table 2. DEGs related to hormone synthesis in response to drought stress.
Table 2. DEGs related to hormone synthesis in response to drought stress.
TermGene IDlog2FCGene DescriptionFDR
ABAPYL/PYRc72499.graph_c25.61abscisic acid receptor PYR/PYL family (A)5.09 × 10−14
c64811.graph_c0−2.31abscisic acid receptor PYR/PYL family (A)2.69 × 10−58
c73702.graph_c1−2.28K14496 abscisic acid receptor PYR/PYL family (A)0.00000743
PP2Cc86830.graph_c02.60probable protein phosphatase 2C 511.62 × 10−19
SAPR1c31398.graph_c04.97basic form of pathogenesis-related protein 1-like2.11 × 10−159
JAJAZc75424.graph_c02.43protein TIFY 10B-like4.90 × 10−58
c75566.graph_c02.14jasmonate ZIM domain-containing protein (A)8.12 × 10−44
c77115.graph_c13.17jasmonate ZIM domain-containing protein (A)2.58 × 10−74
c77115.graph_c23.22Protein TIFY 10B2.80 × 10−87
c88229.graph_c02.13protein TIFY 9-like0.000018
MYC2c88848.graph_c12.14transcription factor MYC2-like1.13 × 10−24
AuxinGH3c78593.graph_c12.24auxin responsive GH3 gene family (A)5.93 × 10−6
c83994.graph_c03.57auxin responsive GH3 gene family (A)3.63 × 10−62
SAURc76579.graph_c02.39uncharacterized protein3.12 × 10−16
c80406.graph_c52.71hypothetical protein MIMGU_mgv1a0212152mg4.51 × 10−26
c63583.graph_c0−4.29auxin-induced protein 10A51.74 × 10−17
c64412.graph_c0−2.18SAUR family protein (A)5.04 × 10−9
c65963.graph_c0−3.66indole-3-acetic acid-induced protein ARG7-like1.28 × 10−13
c84555.graph_c1−4.21SAUR family protein (A)4.31 × 10−11
Ethyle-neMPK6c75482.graph_c02.19mitogen-activated protein kinase 81.2624 × 10−3
EBF1/2c70061.graph_c02.28EIN3-binding F-box protein (A)4.3002 × 10−3
Table 3. DEGs in “Photosynthesis-antenna Proteins” pathway of KEGG.
Table 3. DEGs in “Photosynthesis-antenna Proteins” pathway of KEGG.
TermGene IDlog2FCGene DescriptionFDR
LHCA1c75167.graph_c0−2.30chlorophyll a-b binding protein 6, chloroplastic1.13 × 10−47
LHCA2c57238.graph_c0−2.30chlorophyll a-b binding protein, chloroplastic1.52 × 10−37
LHCA3c71085.graph_c0−2.07chlorophyll a-b binding protein 8, chloroplastic-like3.11 × 10−14
c85515.graph_c0−2.33chlorophyll a-b binding protein 8, chloroplastic2.41 × 10−138
c71085.graph_c1−2.20chlorophyll a-b binding protein 8, chloroplastic-like5.63 × 10−5
LHCA4c57961.graph_c0−3.78chlorophyll a-b binding protein 4, chloroplastic3.42 × 10−11
c81195.graph_c1−3.42agamous-like MADS-box protein AGL21 isoform X31.69 × 10−77
c31746.graph_c0−3.92chlorophyll a-b binding protein P4, chloroplastic-like1.20 × 10−120
LHCB1c85665.graph_c1−3.25chlorophyll a/b-binding protein PS II-Type I6.78 × 10−29
c85665.graph_c2−3.68chlorophyll a-b binding protein 21, chloroplastic-like6.40 × 10−63
c83506.graph_c0−3.53chlorophyll a/b-binding protein, partial1.02 × 10−12
LHCB2c31726.graph_c0−2.58chlorophyll a-b binding protein 5, chloroplastic1.26 × 10−46
c77073.graph_c0−3.51chlorophyll A/B binding protein, putative7.67 × 10−58
LHCB3c82382.graph_c0−2.60chlorophyll a-b binding protein 13, chloroplastic5.22 × 10−47
c82382.graph_c1−2.74chlorophyll a-b binding protein 13, chloroplastic5.88 × 10−32
c84778.graph_c0−2.35chlorophyll a-b binding protein 13, chloroplastic1.14 × 10−20
LHCB4c57394.graph_c0−3.66chlorophyll a-b binding protein CP29.1, chloroplastic4.63 × 10−60
LHCB5c72073.graph_c1−2.03chlorophyll a-b binding protein CP26, chloroplastic4.96 × 10−43
c72073.graph_c0−2.25chlorophyll a-b binding protein CP26, chloroplastic7.71 × 10−40
LHCB6c76630.graph_c1−2.38hypothetical protein MIMGU_mgv1a012260mg2.07 × 10−23
Table 4. DEGs in “Betalain biosynthesis” pathway of KEGG.
Table 4. DEGs in “Betalain biosynthesis” pathway of KEGG.
TermGene IDlog2FCGene DescriptionFDR
TYRc83086.graph_c02.69Tyrosinase2.09 × 105
COMTc26366.graph_c03.02catechol O-methyltransferase7.64 × 105
DOPAc75132.graph_c03.38PREDICTED: 4,5-DOPA dioxygenase extradiol-like1.87 × 1089
Table 5. DEGs in “Flavone and flavonol biosynthese” pathway of KEGG.
Table 5. DEGs in “Flavone and flavonol biosynthese” pathway of KEGG.
TermGene IDlog2FCGene DescriptionFDR
E1.14.13.21c32062.graph_c02.11benzoate 4-monooxygenase cytochrome P4501.06 × 105
AOMTc57467.graph_c04.20PREDICTED: flavonoid 3&apos; 5&apos; -methyltransferase-like1.59 × 1013
c67675.graph_c03.17PREDICTED: flavonoid 3&apos; 5&apos; -methyltransferase-like2.84 × 1022
C12RT1c69454.graph_c04.05hypothetical protein MIMGU_mgv1a022315mg 1.81 × 1017
Table 6. DEGs in “Nitrogen metabolism” pathway of KEGG.
Table 6. DEGs in “Nitrogen metabolism” pathway of KEGG.
TermGene IDlog2FCGene DescriptionFDR
NRc88329.graph_c0−2.65Nitrate reductase 21.32 × 10−97
NirAc85021.graph_c0−2.64Ferredoxin–nitrite reductase4.29 × 10−84
GLULc89561.graph_c02.21glutamine synthetase42.24 × 107
GLT1c85092.graph_c22.29glutamate synthase (NADPH/NADH)1.53 × 105
Table 7. DEGs in “Ubiquitin mediated proteolysis” pathway of KEGG.
Table 7. DEGs in “Ubiquitin mediated proteolysis” pathway of KEGG.
TermGene IDlog2FCGene DescriptionFDR
E2UBE2Ac56569.graph_c02.73ubiquitin-conjugating enzyme E2 A9.53 × 10−7
UBE2Oc78080.graph_c12.32ubiquitin-conjugating enzyme E2 O;
A orthologs to drought gene GmMYB177
1.03 × 10−6
UBE2Wc43734.graph_c03.58ubiquitin-conjugating enzyme E2 W3.64 × 10−6
UBE2Nc61661.graph_c02.30ubiquitin-conjugating enzyme E2 N3.35 × 10−6
UBE2D-Ec89609.graph_c02.54ubiquitin-conjugating enzyme E2 D/E4.15 × 10−6
UBE2Ic46599.graph_c02.71ubiquitin-conjugating enzyme E2 I;
A orthologs to drought gene GmMYB177
3.46 × 10−8
UBE2G1c26287.graph_c02.14ubiquitin-conjugating enzyme E2 G12.449 × 10−3
c25902.graph_c03.11ubiquitin-conjugating enzyme E2 G11.23 × 10−12
E3ARF-BP1c84837.graph_c12.40E3 ubiquitin-protein ligase HUWE12.58614 × 10−4
UBE4Bc71025.graph_c02.32ubiquitin conjugation factor E4 B1.45681 × 10−4
CYC4c75600.graph_c02.46peptidyl-prolyl cis-trans isomerase-like 21.46 × 10−5
PRP19c60805.graph_c02.06pre-mRNA-processing factor 192.56 × 10−5
Cul3c63726.graph_c02.00cullin 3 (A)2.73 × 10−5
CYC4c75600.graph_c02.46peptidyl-prolyl cis-trans isomerase-like 21.46 × 10−5
SYVNc79541.graph_c02.01ubiquitin-protein ligase synoviolin1.091839 × 10−3
Cdh1c47817.graph_c02.93cell division cycle 20-like protein 1, cofactor of APC complex (A)1.53 × 10−5
TRIP12c82561.graph_c02.52E3 ubiquitin-protein ligase TRIP122.08 × 10−9

Share and Cite

MDPI and ACS Style

Wang, B.; Lv, X.-Q.; He, L.; Zhao, Q.; Xu, M.-S.; Zhang, L.; Jia, Y.; Zhang, F.; Liu, F.-L.; Liu, Q.-L. Whole-Transcriptome Sequence Analysis of Verbena bonariensis in Response to Drought Stress. Int. J. Mol. Sci. 2018, 19, 1751. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms19061751

AMA Style

Wang B, Lv X-Q, He L, Zhao Q, Xu M-S, Zhang L, Jia Y, Zhang F, Liu F-L, Liu Q-L. Whole-Transcriptome Sequence Analysis of Verbena bonariensis in Response to Drought Stress. International Journal of Molecular Sciences. 2018; 19(6):1751. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms19061751

Chicago/Turabian Style

Wang, Bei, Xue-Qi Lv, Ling He, Qian Zhao, Mao-Sheng Xu, Lei Zhang, Yin Jia, Fan Zhang, Feng-Luan Liu, and Qing-Lin Liu. 2018. "Whole-Transcriptome Sequence Analysis of Verbena bonariensis in Response to Drought Stress" International Journal of Molecular Sciences 19, no. 6: 1751. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms19061751

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop