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

**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 **Figure 4.** databases and the ordinate is the number of DEGs annotated in it. 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 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. **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.

#### In this study, a total of 2146 (44.44%) transcription factor DEGs were identified, of which 1656 *2.4. DEGs of Transcription Factors (TFs) under Drought Stress 2.4. DEGs of Transcription Factors (TFs) under Drought Stress*

drought stress.

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

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

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

#### *2.5. Expression Level of DEGs' Changes and Verification Using qRT-PCR 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). 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).

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

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

**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−ΔΔ*C*<sup>t</sup> 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−∆∆*C*<sup>t</sup> method.

#### **3. Discussion 3. Discussion**

#### *3.1. Morphological and Physiological Index Analysis 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 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.

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