*2.6. Expression Patters of EgrVQ Genes in Different Tissues of E. grandis*

In our study, the expression patterns of *EgrVQ* genes were examined in different tissues, including root, xylem, phloem, mature leaves, and young leaves (Figure 6). The majority (88.9%) of *EgrVQ* genes were expressed in all tissues. However, almost all *EgrVQ*s showed differential gene expression. For example, from subfamily I members, *EgrVQ2* and *EgrVQ22* were relative highly expressed when compared to other members, *EgrVQ13* and *EgrVQ14* mainly were expressed in leaves, and *EgrVQ18* was not expressed in the root. Moreover, from subfamily V members, when compared to *EgrVQ20* and *EgrVQ21, EgrVQ9* and *EgrVQ11* were relative highly expressed in most tissues. In addition, from subfamily VII members, *EgrVQ1* and *EgrVQ26* were very highly expressed in root and xylem tissues than other tissues. Therefore, *EgrVQ1* and *EgrVQ26* might participate in root development and play vital roles in the xylem development of *E. grandis*.

When considering that 92.6% and 88.9% of *EgrVQ* genes were expressed in young leaves and mature leaves, respectively, we selected leaves for determining the expression profiles of *EgrVQ* genes in response to different treatments. *Int. J. Mol. Sci.* **2019**, *20*, x 9 of 17

**Figure 6.** Semi-quantitative Real-Time PCR (RT-PCR) of the expression of *EgrVQ* genes in different tissues. From left to right: root, xylem, phloem, mature leaves, and young leaves. All of these tissues were collected from 6-week-old GL1 plants. **Figure 6.** Semi-quantitative Real-Time PCR (RT-PCR) of the expression of *EgrVQ* genes in different tissues. From left to right: root, xylem, phloem, mature leaves, and young leaves. All of these tissues were collected from 6-week-old GL1 plants.

#### *2.7. Expression Profiles of EgrVQ Genes in Response to Hormones*

up-regulated.

treatment.

*2.7. Expression Profiles of EgrVQ Genes in Response to Hormones*  To investigate the function of *EgrVQ* genes in response to biotic stresses, we determined the expression patterns of *EgrVQ* genes under various plant hormone treatments, including BRs, MeJA, SA, and ABA (Figure 7). The results showed that the *EgrVQ* genes were expressed in diverse patterns under different hormone treatments. Most of the *EgrVQ* genes presented differential expression under the different phytohormone treatments. Some differences are extremely To investigate the function of *EgrVQ* genes in response to biotic stresses, we determined the expression patterns of *EgrVQ* genes under various plant hormone treatments, including BRs, MeJA, SA, and ABA (Figure 7). The results showed that the *EgrVQ* genes were expressed in diverse patterns under different hormone treatments. Most of the *EgrVQ* genes presented differential expression under the different phytohormone treatments. Some differences are extremely significant when compared to the reference gene *EgrEF*.

significant when compared to the reference gene *EgrEF*. In the BR treatment, most of *EgrVQ* genes were up-regulated and the *EgrVQs* were clustered into three main groups according to their expression patterns (Figure 7A). Group BR-1 contained three *EgrVQs,* including *EgrVQ19*, *EgrVQ21*, and *EgrVQ27*, which were significantly decreased with the increasing duration of BR treatment. There were three *EgrVQs* in group BR-2 and their expression levels were dramatically increased at 1 h and 168 h after BR treatment.The rest of the In the BR treatment, most of *EgrVQ* genes were up-regulated and the *EgrVQs* were clustered into three main groups according to their expression patterns (Figure 7A). Group BR-1 contained three *EgrVQs*, including *EgrVQ19*, *EgrVQ21*, and *EgrVQ27*, which were significantly decreased with the increasing duration of BR treatment. There were three *EgrVQs* in group BR-2 and their expression levels were dramatically increased at 1 h and 168 h after BR treatment.The rest of the *EgrVQ* genes were included in group BR-3, which contained 19 *EgrVQs* and they were continuously up-regulated.

*EgrVQ* genes were included in group BR-3, which contained 19 *EgrVQs* and they were continuously

significantly altered at 1 h to 24 h after SA treatment. This result indicated that the expression of most *EgrVQs* was affected by long-term treatment with SA. The rest of the *EgrVQ*s were clustered in

group SA-2 and their expression was significantly reduced or showed no obvious change.

Under SA treatment, the *EgrVQs* were mainly clustered into two groups according to their expression patterns (Figure 7B). The majority of *EgrVQs* (22/27) were clustered into group SA-1, in

In the MeJA treatment, the *EgrVQs* were mainly clustered into three groups according to their expression patterns (Figure 7C). In the MeJA-1 group, there were seven *EgrVQs,* and their expressions were greatly decreased at 1 h, 24 h, and 168 h after MeJA treatment. The MeJA-2 group contained eleven *EgrVQs*, which showed significantly decreased expression at 1 h and then increased the expression at 6 h and reached their maximum expression level and decreased afterward. The rest of *EgrVQ*s were clustered in group MeJA-3, in which the expression levels of *EgrVQs* were greatly down-regulated at 1–24 h, and increased back to the level of the control or increased. However, it was found that *EgrVQ*s were immediately down-regulated under MeJA

Under SA treatment, the *EgrVQs* were mainly clustered into two groups according to their expression patterns (Figure 7B). The majority of *EgrVQs* (22/27) were clustered into group SA-1, in which the expression levels of *EgrVQs* were greatly up-regulated at 168 h, but they were not significantly altered at 1 h to 24 h after SA treatment. This result indicated that the expression of most *EgrVQs* was affected by long-term treatment with SA. The rest of the *EgrVQ*s were clustered in group SA-2 and their expression was significantly reduced or showed no obvious change. Under ABA treatment, the expression level of the majority of *EgrVQs* was obviously reduced at 168 h and the *EgrVQs* were mainly clustered into two groups (Figure 7D). In group ABA-1, the expression of *EgrVQs* was significantly up-regulated at 1 h and then decreased with increasing treatment time under ABA treatment. In group ABA-2, there were nine *EgrVQs*, and their expression levels increased at 1 h, 6 h, or 24 h.

*Int. J. Mol. Sci.* **2019**, *20*, x 10 of 17

**Figure 7.** Expression analysis of 25 *EgrVQ* genes in *E. grandis* following brassinosteroid (BR) (**A**), salicylic acid (SA) (**B**), methyl jasmonate (MeJA) (**C**), and abscisic acid (ABA) (**D**) treatment, as determined by Real-Time quantitative PCR (qRT-PCR). The relative expression levels were calculated using the 2−ΔΔCt method. The heatmap was created using MEV. **Figure 7.** Expression analysis of 25 *EgrVQ* genes in *E. grandis* following brassinosteroid (BR) (**A**), salicylic acid (SA) (**B**), methyl jasmonate (MeJA) (**C**), and abscisic acid (ABA) (**D**) treatment, as determined by Real-Time quantitative PCR (qRT-PCR). The relative expression levels were calculated using the 2−∆∆Ct method. The heatmap was created using MEV.

*2.8. Expression Profiles of EgrVQ Genes in Response to Stress Treatments*  To research the function of *EgrVQ* genes in response to abiotic stressors, we investigated the expression patterns of *EgrVQ* genes under NaCl, cold, and heat treatments, respectively (Figure 8). Overall, the *EgrVQ* genes showed differential expression under the different abiotic stresses. Some differences are even significant when compared to the reference gene *EgrEF*. Specifically, most of In the MeJA treatment, the *EgrVQs* were mainly clustered into three groups according to their expression patterns (Figure 7C). In the MeJA-1 group, there were seven *EgrVQs*, and their expressions were greatly decreased at 1 h, 24 h, and 168 h after MeJA treatment. The MeJA-2 group contained eleven *EgrVQs*, which showed significantly decreased expression at 1 h and then increased the expression at 6 h and reached their maximum expression level and decreased afterward. The rest of *EgrVQ*s were clustered in group MeJA-3, in which the expression levels of *EgrVQs* were greatly down-regulated at

*EgrVQ* genes were up-regulated at 1 h or 6 h and they reached their highest expression levels, but they were then down-regulated at 168 h. These results indicated that *EgrVQ* genes could respond to short-term temperature stress. Interestingly, it was found that similar expression patterns of *EgrVQ*  1–24 h, and increased back to the level of the control or increased. However, it was found that *EgrVQ*s were immediately down-regulated under MeJA treatment.

Under ABA treatment, the expression level of the majority of *EgrVQs* was obviously reduced at 168 h and the *EgrVQs* were mainly clustered into two groups (Figure 7D). In group ABA-1, the expression of *EgrVQs* was significantly up-regulated at 1 h and then decreased with increasing treatment time under ABA treatment. In group ABA-2, there were nine *EgrVQs*, and their expression levels increased at 1 h, 6 h, or 24 h. *Int. J. Mol. Sci.* **2019**, *20*, x 11 of 17 genes occurred under cold and heat treatments.Subsequently, *EgrVQs* were classified into three groups both in the cold and heat treatments, respectively. In group Cold-1 and group Heat-2, there were 10 and 15 *EgrVQs*, respectively. The expression

#### *2.8. Expression Profiles of EgrVQ Genes in Response to Stress Treatments* patterns of these genes were similar, and their expressions were greatly increased at 6-48 h and then decreased after ABA treatment. Similarly, 10 and seven *EgrVQs* were classified into group Cold-2

To research the function of *EgrVQ* genes in response to abiotic stressors, we investigated the expression patterns of *EgrVQ* genes under NaCl, cold, and heat treatments, respectively (Figure 8). Overall, the *EgrVQ* genes showed differential expression under the different abiotic stresses. Some differences are even significant when compared to the reference gene *EgrEF*. Specifically, most of *EgrVQ* genes were up-regulated at 1 h or 6 h and they reached their highest expression levels, but they were then down-regulated at 168 h. These results indicated that *EgrVQ* genes could respond to short-term temperature stress. Interestingly, it was found that similar expression patterns of *EgrVQ* genes occurred under cold and heat treatments.Subsequently, *EgrVQs* were classified into three groups both in the cold and heat treatments, respectively. and group Heat-1, respectively. Their expressions levels were greatly increased at 1 h and 6 h, but they then decreased after ABA treatment. In group Cold-3, *EgrVQs* were significantly increased at 1 h. In group Heat-3, *EgrVQs* were significantly increased at 48 h. Under NaCl treatment, the *EgrVQs* were clustered into three groups according to their expression patterns (Figure 8C). In group NaCl-1, the expression of *EgrVQs* showed slight decreases with an increasing treatment time. In group NaCl-2, *EgrVQ*s were down-regulated at 1 h and then up-regulated with increasing NaCl treatment duration at 6 h and 24 h. The remaining of *EgrVQs* were classified to group NaCl-3. In this group, the expression of *EgrVQs* gradually increased and reached their highest expression level at 168 h.

**Figure 8.** Expression analysis of 25 *EgrVQ* genes in *E. grandis* following cold (**A**), heat (**B**), and NaCl (**C**) treatments, as determined by qRT-PCR. The relative expression levels were calculated using the **Figure 8.** Expression analysis of 25 *EgrVQ* genes in *E. grandis* following cold (**A**), heat (**B**), and NaCl (**C**) treatments, as determined by qRT-PCR. The relative expression levels were calculated using the 2 <sup>−</sup>∆∆Ct method. The heatmap was created using MEV.

**3. Discussion**  In group Cold-1 and group Heat-2, there were 10 and 15 *EgrVQs*, respectively. The expression patterns of these genes were similar, and their expressions were greatly increased at 6-48 h and then

In our study, 27 *EgrVQ* genes were identified using the genome database in *E. grandis*. Subsequently, phylogenetic analysis, conserved motifs, and analysis of *cis*-elements of *EgrVQs* were

2−ΔΔCt method. The heatmap was created using MEV.

decreased after ABA treatment. Similarly, 10 and seven *EgrVQs* were classified into group Cold-2 and group Heat-1, respectively. Their expressions levels were greatly increased at 1 h and 6 h, but they then decreased after ABA treatment. In group Cold-3, *EgrVQs* were significantly increased at 1 h. In group Heat-3, *EgrVQs* were significantly increased at 48 h.

Under NaCl treatment, the *EgrVQs* were clustered into three groups according to their expression patterns (Figure 8C). In group NaCl-1, the expression of *EgrVQs* showed slight decreases with an increasing treatment time. In group NaCl-2, *EgrVQ*s were down-regulated at 1 h and then up-regulated with increasing NaCl treatment duration at 6 h and 24 h. The remaining of *EgrVQs* were classified to group NaCl-3. In this group, the expression of *EgrVQs* gradually increased and reached their highest expression level at 168 h.

#### **3. Discussion**

In our study, 27 *EgrVQ* genes were identified using the genome database in *E. grandis*. Subsequently, phylogenetic analysis, conserved motifs, and analysis of *cis*-elements of *EgrVQs* were performed. All of the *EgrVQs* encoded relatively small proteins of less than 400 amino acids, and most VQ motif-containing proteins contained FxxVQxLTG(20/27), FxxVQxFTG(4/27), and FxxVQxVTG(2/27). Moreover, most *EgrVQs*(25/27) only contained an exon. Our results were consistent with *A. thaliana*, poplar, and rice [7,17,20], which indicated that *VQ* genes were relatively conservative. However, motif FxxVQxLSG(1/27), which was not shown in *A. thaliana*, rice, poplar, or grape, was found in *E. grandis* (Table 3).

Additionally, according to 20 conserved motifs in the analysis of *EgrVQ* genes with Multiple Expectation Maximization for Motif Elicitation (MEME), motif1 was presented in all of the *EgrVQ* genes (Figure 5). This result was consistent with other plants, like poplar [7] and moso bamboo [8,27], implying that the motif 1 imparts specific functions to the VQ protein. Moreover, it was found that several *cis*-elements that related to plant hormones and abiotic stressors were identified in the promoter region of the *EgrVQ* genes (Table 4). Similarly, most of these *cis*-elements occurred in other plants, like rice [17], poplar [7], and bamboo [8,27]. Interestingly, DRE is a *cis*-acting element that is involved in salt stress,cold stress, and dehydration, which was only presented in *EgrVQ20* and *EgrVQ21* (Table 4), and also seldom occurred in other plants.

In a previous study, VQ proteins participated in regulating diverse developmental processes, especially in response to biotic and abiotic stressors [20]. SA and JA are vital defense signaling molecules to response to pathogen infection and other abiotic stress conditions. In our study, the transcription of most *EgrVQ* genes was altered under SA and MeJA treatments. Interestingly, the expression of many *EgrVQ* genes was significantly increased at 168 h of SA treatment, which was not studied in moso bamboo [27], *A. thaliana* [20], and *Vitisvinifera* [28]. It is worth noting that *EgrVQ2*, *EgrVQ18*, and *EgrVQ22*, which were homologs of *AtVQ23* and *AtVQ16*, were highly expressed at 168 h under SA treatment. In *A. thaliana*, the over-expression of *AtVQ23* induced hyper-activate defense-related genes in plants following pathogen infection or SA and MeJA treatments, which enhanced resistance to infection by *Pseudomonassyringae* [29]. Meanwhile, *AtVQ16* could regulate the immune response by regulating *AtWRKY33* and further stimulating the DNA binding activity [30]. It was also found that *EgrVQ2*2 expression also increased under MeJA treatment. However, *EgrVQ2* expression was not significantly changed and *EgrVQ18* expression was slightly down-regulated under MeJA treatment. Hence, these results implied that *EgrVQ18* could respond to SA treatment and it was further involved in SA-mediated defense responses by the long-term effect. *EgrVQ22* could also respond to SA and MeJA treatments.

Moreover, *AtVQ22*, which was positively regulated by the COI1 (CORONATINEINSENSITIVE1) dependent signaling pathway, was a master controller regulating JA-mediated plant response against insects and pathogens [16,31]. *EgrVQ8*, which was homologous with *AtVQ22*, was induced at 6 h under MeJA treatment and 168 h under SA treatment, indicating that *EgrVQ8* has a similar function in

*E. grandis.* However, *EgrVQ26*, which is another homologous gene with *AtVQ22*, was not changed under MeJA treatment.

Furthermore, ABA and BR were also important in the response of plants to abiotic stressors [32,33]. In our study, most *EgrVQs* were up-regulated at 1 h–24 h when the plants were exposed to ABA, which was similar to moso bamboo *VQ* genes [8]. Meanwhile, only three genes showed up-regulation and multiple *OsVQ* genes showed down-regulation after 12 h under ABA treatment [17]. However, we found that 22 *EgrVQs* showed significantly down-regulation at 168 h under ABA treatment. These results indicated that the VQ genes had special roles in response to ABA treatment in *E. grandis*. It was also found that *EgrVQ* subfamily I members, including *EgrVQ2*, *EgrVQ13*, *EgrVQ14*, *EgrVQ22*,and *EgrVQ18*, were highly expressed under ABA treatment,which was similar with *OsVQ6*, *OsVQ7*, and *OsVQ10* belonging to subfamily. These results implied that these genes have similar functions and they might be related to the ABA signaling pathway. Specifically, *AtVQ15* regulated the osmotic stress tolerance with *WRKY25* and *WRKY51* in *A. thaliana* [20]. Its homologous gene, *EgrVQ27*, was significantly up-regulated after the ABA treatment, suggesting that *EgrVQ27* might affect osmotic stress tolerance in *E. grandis*. In addition, we found that most of the *EgrVQ* genes were highly expressed under the BR treatment (Figure 7A). BRASSINAZOLE-RESISTANT 1(BZR1), which is a key transcription factor (TF) in the BRs signal pathway, can trigger the expression of *FLS2* and *SNC1* to enhance pathogen resistance in plants [34,35].These results implied that *EgrVQ* genes might take part in enhancing the pathogen resistance of the BR signal pathway.

Plants also need to evolve efficient defense systems to protect themselves from abiotic stressors, like extreme temperature, salt, and so on. AtVQ9 acted as a repressor of WRKY8 to maintain a balance in the activation of WRKY8-mediated signaling pathways that are involved in salt stress in *A. thaliana* [22]. In our study, *EgrVQ5* and *EgrVQ9*, both homologous of *AtVQ9*, showed the same expression trends and they were increased at 168 h under NaCl treatment (Figure 8C). This indicated that the *VQ* genes were relatively conservative and were also involved in salt stress. In contrast, another homologous gene, *EgrVQ11*, was decreased at 168 h. *AtVQ15* was induced by dehydration and high salinity, whereas lines that overexpressed this gene showed an increased desensitivity to salt stress during seed germination and seedling growth [20,22]. Interestingly, its homologous gene, *EgrVQ1*, showed obviously decreased expression under NaCl treatment. We speculate that there is a difference in the woody and herbaceous plant.

Furthermore, few previous studies have been conducted on the function of *VQ* genes under heat and cold stresses in the woody plants. Our results showed that 23 and 24 *EgrVQ* genes were up-regulated at 1 h or 6 h, respectively, and then down-regulated at 168 h under cold and heat treatments, respectively (Figure 8A,B), which was consistent with Chinese Cabbage (*Brassica pekinensis*) [36]. It is noteworthy that most of the *EgrVQ* genes showed similar expression patterns between heat and cold treatment, suggesting that *EgrVQ*s existed in the same response pathway to heat and cold stress.

Overall, our results implied various stressors motivated or repressed a batch of *EgrVQ*s. These results might aid in the selection of available candidate genes in the *EgrVQ* gene family for deeper functional characterization.

#### **4. Materials and Methods**

#### *4.1. Identification of VQ Gene Family in E. grandis*

We conducted a systematic search on related bioinformatics analysis websites in order to acquire all of the proteins in *Eucalyptus*. Firstly, the sequences of *VQ* genes of *A. thaliana*, rice, and poplar were downloaded. Afterwards, the *VQ* genes were searched in *E. grandis* was using in NCBI (https:// www.ncbi.nlm.nih.gov/), UniProt (https://www.uniprot.org), and EucGenIE (https://eucgenie.org/) databases. Next, all of the predicted *VQ* genes were subjected to SMART (http://smart.emblheidelberg. de/) and Pfam (https://pfam.sanger.ac.uk/) to confirm that they contained the VQ motif (PF05678). Lastly, physical parameters of *VQ* genes, including open reading frame (ORF) length, protein length, isoelectric point (pI), and molecular weight were calculated in ExPASy (http://www.expasy.org/ tools) [37]. Subcellular localization was predicted using the WoLF PSORT (http://wolfpsort.org/) [38] and TargetP 1.1 (http://www.cbs.dtu.dk/services/TargetP/) servers [39].

## *4.2. Chromosomal Location, Gene Duplication, and Identification of Paralogs and Orthologs of EgrVQ Genes*

MapChart 2.30 software drew the picture of chromosomal location [40] on the basis of initial position information that was provided in Phytozomev12.1.6 (https://phytozome.jgi.doe.gov/pz/ portal.html#). In addition, *EgrVQ* genes, which are located on duplicated chromosomal blocks, were considered to undergo segmental duplication [26]. Next, we used a previously described method for the analysis of paralogs and orthologs [41]. For the elected species, all-against-all nucleotide sequence similarity searches were done among the protein sequences by BLASTN software [42]. Specifically, in one species, a pair of matching sequences was shown as pairs of paralogs, which, when aligned, exceed 300bp and the identity was over 40% [7]. In two different species, the two sequences were defined as orthologs with the reciprocal best hits for each being within >300bp of the aligned sequences from two different species [7,8].
