**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–9], *Oryza sativa* L. [10–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. 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–9], *Oryza sativa* L. [10–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. Results**

#### *2.1. Phenotypic and Physiological Indicators of Verbena under Drought Stress 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 (Figures 2d and 3a,b), malonaldehyde (MDA) content also increased and relative water content (RWC) did not significantly decrease (Figure 3c,d). 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 (Figures 2d and 3a,b), malonaldehyde (MDA) content also increased and relative water content (RWC) did not significantly decrease (Figure 3c,d).

**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.

*Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 3 of 20

*Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 3 of 20

**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. **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 3.** *Cont.*

*Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 4 of 20

**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.

#### *2.2. Sequencing and Annotation of Transcription and Unigenes 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). 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 × 10−<sup>5</sup> and HMMER parameters E-value ≤ 1 × 10−<sup>10</sup> 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. 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 <sup>≤</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> and HMMER parameters E-value <sup>≤</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>10</sup> 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). 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), 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).

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 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 higherorder 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).
