**2. Results**

#### *2.1. Identification and Sequence Analysis of the PtrMORF Gene Family* to drought stress in poplar.

We searched the poplar genome with known *A. thaliana* MORF proteins as queries. Initially, nine putative *PtrMORF* genes were obtained—*PtrMORF1.1*, *PtrMORF1.2*, *PtrMORF1.3*, *PtrMORF2.1*, *PtrMORF2.2*, *PtrMORF3*, *PtrMORF8.1*, *PtrMORF8.2*, and *PtrMORF9* based on a phylogenetic analysis using the amino acid sequences of all MORFs of *P. trichocarpa* as well as those of *A. thaliana* (Figure 1A). Poplar MORF proteins were predicted using TargetP and Wolf PSORT to enter mitochondria, chloroplasts, or nuclei, like their homologs in *A. thaliana* (Figure 1A). No known motif in poplar MORF proteins was found in the PFAM and INTERPRO databases, but the MORF box was identified, as in previous studies [27,28,33]. Novel putative motifs were explored using the MEME server. By selecting a motif length of between 15 and 50 aa, we identified four conserved motifs located in the MORF domain (Figure 1B–D). **2. Results**  *2.1. Identification and Sequence Analysis of the PtrMORF Gene Family*  We searched the poplar genome with known *A. thaliana* MORF proteins as queries. Initially, nine putative *PtrMORF* genes were obtained—*PtrMORF1.1*, *PtrMORF1.2*, *PtrMORF1.3*, *PtrMORF2.1*, *PtrMORF2.2*, *PtrMORF3*, *PtrMORF8.1*, *PtrMORF8.2*, and *PtrMORF9*based on a phylogenetic analysis using the amino acid sequences of all MORFs of *P. trichocarpa* as well as those of *A. thaliana* (Figure 1A). Poplar MORF proteins were predicted using TargetP and Wolf PSORT to enter mitochondria, chloroplasts, or nuclei, like their homologs in *A. thaliana* (Figure 1A). No known motif in poplar MORF proteins was found in the PFAM and INTERPRO databases, but the MORF box was identified, as in previous studies [27,28,33]. Novel putative motifs were explored using the MEME server. By selecting a motif length of between 15 and 50 aa, we identified four conserved motifs located in the MORF domain (Figure 1B–D).

**Figure 1.** Poplar multiple organellar RNA editing factor (MORF) proteins and their conserved motifs. (**A**) Multiple sequences alignment of *A. thaliana* and *P. trichocarpa* MORF proteins was carried out using MUSCLE, and the Neighbor-Joining (NJ) tree was built using MEGA v7.0. And the chloroplast, mitochondrial or nuclei transit peptides of poplar MORF proteins were predicted using TargetP (http://www.cbs.dtu.dk/services/TargetP/) and Wolf PSORT (http://wolfpsort.org/). M, mitochondria; C, chloroplast; Nucl, nuclei. Green diamond, MORF9 and PtrMORF9; Blue triangle, MORF1, PtrMORF1.1, PtrMORF1.2 and PtrMORF1.3; Red circle, MORF2, PtrMORF2.1 and PtrMORF2.2; **Figure 1.** Poplar multiple organellar RNA editing factor (MORF) proteins and their conserved motifs. (**A**) Multiple sequences alignment of *A. thaliana* and *P. trichocarpa* MORF proteins was carried out using MUSCLE, and the Neighbor-Joining (NJ) tree was built using MEGA v7.0. And the chloroplast, mitochondrial or nuclei transit peptides of poplar MORF proteins were predicted using TargetP (http://www.cbs.dtu.dk/services/TargetP/) and Wolf PSORT (http://wolfpsort.org/). M, mitochondria; C, chloroplast; Nucl, nuclei. Green diamond, MORF9 and PtrMORF9; Blue triangle, MORF1, PtrMORF1.1, PtrMORF1.2 and PtrMORF1.3; Red circle, MORF2, PtrMORF2.1 and PtrMORF2.2; Orange circle, MORF3 and PtrMORF3; Orange diamond, MORF8, PtrMORF8.1 and PtrMORF8.2. Local meant the gene where was located in (**B**). Alignment of conserved MORF domains in poplar MORF proteins was conducted using DNAMAN (https://www.lynnon.com/). Lake blue, Motif1; purple, Motif2; red, Motif3; green, Motif4. (**C**,**D**) Putative motifs were explored using the MEME server with the parameters of between 15 and 50 aa in length and sharing of each motif among all PtrMORF proteins.

We further analyzed the sequence structures of the nine *PtrMORF* genes. An alignment of the genomic sequences to predicted CDS sequences of *PtrMORF* genes, showed that *PtrMORF* genes had a conserved gene structure. The *PtrMORF* genes had three introns with intron phases 0, 1, and 2. In all *PtrMORF* genes, motif 1 was encoded by exon 1, motif 2 was encoded by exon1 and exon 2, and motif 3 was encoded by exons 3 and 4, but motif 4 was only located in exon 4 of *PtrMORF1.1* (Figure 2A,B). The nine *PtrMORF* genes ranged from 1515 to 4441 bp and contained four or five exons (Figure 2A). All identified poplar *MORF* genes encoded proteins ranging from 229 to 470 amino acids, and their sequences contained zero to two transmembrane domains (TMDs). The molecular weight (MW) of the nine putative proteins ranged from 26.0 to 51.7 kDa. The GRAVY values of putative MORFs were negative and ranged from −1.382 to −0.602 (Table S1). We further analyzed the sequence structures of the nine *PtrMORF* genes. An alignment of the genomic sequences to predicted CDS sequences of *PtrMORF* genes, showed that *PtrMORF* genes had a conserved gene structure. The *PtrMORF* genes had three introns with intron phases 0, 1, and 2. In all *PtrMORF* genes, motif 1 was encoded by exon 1, motif 2 was encoded by exon1 and exon 2, and motif 3 was encoded by exons 3 and 4, but motif 4 was only located in exon 4 of *PtrMORF1.1* (Figure 2A,B). The nine *PtrMORF* genes ranged from 1515 to 4441 bp and contained four or five exons (Figure 2A). All identified poplar *MORF* genes encoded proteins ranging from 229 to 470 amino acids, and their sequences contained zero to two transmembrane domains (TMDs). The molecular weight (MW) of the nine putative proteins ranged from 26.0 to 51.7 kDa. The GRAVY values of putative MORFs were negative and ranged from −1.382 to −0.602 (Table S1).

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

Orange circle, MORF3 and PtrMORF3; Orange diamond, MORF8, PtrMORF8.1 and PtrMORF8.2. Local meant the gene where was located in (**B**). Alignment of conserved MORF domains in poplar MORF proteins was conducted using DNAMAN (https://www.lynnon.com/). Lake blue, Motif1; purple, Motif2; red, Motif3; green, Motif4. (**C**,**D**) Putative motifs were explored using the MEME server with the parameters of between 15 and 50 aa in length and sharing of each motif among all

**Figure 2.** Gene structures of poplar *MORF* genes. (**A**,**B**) The gene structures of *PtrMORF* genes were built using GSDraw (http://wheat.pw.usda.gov/piece/GSDraw.php) by submitting both genomic sequences obtained from PopGenIE (http://popgenie.org/) and coding sequences (CDS) of *PtrMORF* **Figure 2.** Gene structures of poplar *MORF* genes. (**A**,**B**) The gene structures of *PtrMORF* genes were built using GSDraw (http://wheat.pw.usda.gov/piece/GSDraw.php) by submitting both genomic sequences obtained from PopGenIE (http://popgenie.org/) and coding sequences (CDS) of *PtrMORF* genes.

#### genes. *2.2. Phylogenetic Comparison of the MORF from Different Species*

*2.2. Phylogenetic Comparison of the MORF from Different Species*  To investigate their molecular evolution and functions in poplar, a phylogenetic analysis of MORF proteins was performed. We used the HMMER 3.0 package to build a Hidden Markov Model (HMM) file (morf.hmm) with 19 MORF domain sequences of MORF proteins from *A. majus*, *P. trichocarpa*, and *A. thaliana* (Table S2). To mine additional MORF domain-encoding genes in other plants, we used the morf.hmm algorithm to query the genomes of six species, representing major evolutionary lineages, including *Arabidopsis lyrata*, *Brachypodium distachyon*, *Glycine max*, *Oryza sativa Japonica*, *Prunus persica*, and *Vitis vinifera*. The numbers of MORF genes in eight species were To investigate their molecular evolution and functions in poplar, a phylogenetic analysis of MORF proteins was performed. We used the HMMER 3.0 package to build a Hidden Markov Model (HMM) file (morf.hmm) with 19 MORF domain sequences of MORF proteins from *A. majus*, *P. trichocarpa*, and *A. thaliana* (Table S2). To mine additional MORF domain-encoding genes in other plants, we used the morf.hmm algorithm to query the genomes of six species, representing major evolutionary lineages, including *Arabidopsis lyrata*, *Brachypodium distachyon*, *Glycine max*, *Oryza sativa Japonica*, *Prunus persica*, and *Vitis vinifera*. The numbers of MORF genes in eight species were comparable, ranging from six (in *P. persica*) to 11 (in *A. lyrata* and *G. max*).

comparable, ranging from six (in *P. persica*) to 11 (in *A. lyrata* and *G. max*). In total, 69 MORF genes were identified in eight plant genomes to build an unrooted tree using MEGA7.0 by employing the neighbor-joining (NJ) method. As shown in Figure 3, sequences were classified into six groups (I, II, III, IV, V, and VI). Each group included *MORF* genes from diverse In total, 69 MORF genes were identified in eight plant genomes to build an unrooted tree using MEGA7.0 by employing the neighbor-joining (NJ) method. As shown in Figure 3, sequences were classified into six groups (I, II, III, IV, V, and VI). Each group included *MORF* genes from diverse plant taxa. Among these, classes III and VI were larger than the others, containing 30 members and accounting for 42.5% of all predicted *MORF* genes. *PtrMORF* genes were found in all classes other than group III.

plant taxa. Among these, classes III and VI were larger than the others, containing 30 members and

**Figure 3.** Phylogenetic relationships of the *MORF* gene family members from poplar (*PtrMORFs*) *Arabidopsis lyrata, Arabidopsis thaliana*, *Brachypodium distachyon*, *Glycine max*, *Oryza sativa Japonica*, *Prunus persica*, and *Vitis vinifera.* The phylogenetic tree was constructed with MEGA 7.0 (http://www.megasofware.net/mega.html/) program by the neighbor-joining method. **Figure 3.** Phylogenetic relationships of the *MORF* gene family members from poplar (*PtrMORFs*) *Arabidopsis lyrata*, *Arabidopsis thaliana*, *Brachypodium distachyon*, *Glycine max*, *Oryza sativa Japonica*, *Prunus persica*, and *Vitis vinifera*. The phylogenetic tree was constructed with MEGA 7.0 (http://www. megasofware.net/mega.html/) program by the neighbor-joining method.

#### *2.3. Chromosomal Distribution, Synteny, and Evolution of PtrMORF Genes 2.3. Chromosomal Distribution, Synteny, and Evolution of PtrMORF Genes*

Figure 4 showed that the genes were distributed on six poplar chromosomes, including chromosomes 1, 3, 4, 8, 10, and 11. Half of the chromosomes had two *PtrMORF* members and the other half had a single *MORF*. Tandem duplication was defined as different members of the gene family occurring within the same or neighboring intergenic region [34]. One tandem duplication event involving two *MORF* genes (*PtrMORF1.2*/*PtrMORF1.3*) was identified by BLASTP and MCScanX methods. In addition to the tandem duplication, seven *PtrMORF* genes were assigned to three segmental duplication events (*PtrMORF1.1/PtrMORF1.2/PtrMORF1.3*, *PtrMORF2.1/PtrMORF2.2*, and *PtrMORF8.1/PtrMORF8.2*) in Populus linkage groups 1, 3, 4, 8, 10, and 11. Notably, *PtrMORF1.2* and *PtrMORF1.3* occurred in both tandem and segmental duplications. A cross-matching event was also found in three doubling blocks of *PtrMORF* genes; the chromosomal fragment in which a gene was located was identical to more than one nonself chromosomal segment. *PtrMORF1.1*, *PtrMORF1.2*, and *PtrMORF1.3* synteny blocks corresponded to the third, fourth, and fourth chromosomes, respectively (Figure 5A and Table S3). Figure 5B showed five of the nine *MORF* genes were involved in three segmental duplication events (*MORF1/MORF4*, *MORF5/MORF6*, and *MORF8*/AT1G53260) in *A. thaliana*. It was worth noting that the gene, AT1G53260, was not included in the scope of *MORF* genes because of partial MORF box in our research [27]. The *MORF* gene of *A. thaliana* was also shown to be highly segmental duplicated, which was similar to that in poplar. Figure 4 showed that the genes were distributed on six poplar chromosomes, including chromosomes 1, 3, 4, 8, 10, and 11. Half of the chromosomes had two *PtrMORF* members and the other half had a single *MORF*. Tandem duplication was defined as different members of the gene family occurring within the same or neighboring intergenic region [34]. One tandem duplication event involving two *MORF* genes (*PtrMORF1.2*/*PtrMORF1.3*) was identified by BLASTP and MCScanX methods. In addition to the tandem duplication, seven *PtrMORF* genes were assigned to three segmental duplication events (*PtrMORF1.1/PtrMORF1.2/PtrMORF1.3*, *PtrMORF2.1/PtrMORF2.2*, and *PtrMORF8.1/PtrMORF8.2*) in Populus linkage groups 1, 3, 4, 8, 10, and 11. Notably, *PtrMORF1.2* and *PtrMORF1.3* occurred in both tandem and segmental duplications. A cross-matching event was also found in three doubling blocks of *PtrMORF* genes; the chromosomal fragment in which a gene was located was identical to more than one nonself chromosomal segment. *PtrMORF1.1*, *PtrMORF1.2*, and *PtrMORF1.3* synteny blocks corresponded to the third, fourth, and fourth chromosomes, respectively (Figure 5A and Table S3). Figure 5B showed five of the nine *MORF* genes were involved in three segmental duplication events (*MORF1/MORF4*, *MORF5/MORF6*, and *MORF8*/AT1G53260) in *A. thaliana*. It was worth noting that the gene, AT1G53260, was not included in the scope of *MORF* genes because of partial MORF box in our research [27]. The *MORF* gene of *A. thaliana* was also shown to be highly segmental duplicated, which was similar to that in poplar.

Mbp in *A. thaliana*.

= 1 supports neutral evolution.

**Figure 4.** Chromosomal locations of Populus *MORF* family members. The chromosomal locations of the *MORF* genes were mapped with MapDraw program. The number to the left of each chromosome represented the size of the chromosome in Mbp. **Figure 4.** Chromosomal locations of Populus *MORF* family members. The chromosomal locations of the *MORF* genes were mapped with MapDraw program. The number to the left of each chromosome represented the size of the chromosome in Mbp. **Figure 4.** Chromosomal locations of Populus *MORF* family members. The chromosomal locations of the *MORF* genes were mapped with MapDraw program. The number to the left of each chromosome represented the size of the chromosome in Mbp.

**Figure 5.** Interchromosomal relationships of *MORF* genes. (**A**) Interchromosomal relationships of *MORF* genes in poplar. (**B**) Interchromosomal relationships of *MORF* genes in *A. thaliana*. Black lines indicated all synteny blocks in the Populus or *A. thaliana* genome and the red indicated synteny blocks where *MORF* duplicated gene pairs were. The chromosome number was indicated at the bottom of each chromosome and each cell on the outside of chromosomes represented 5 Mbp in *Populus* and 1.5 **Figure 5.** Interchromosomal relationships of *MORF* genes. (**A**) Interchromosomal relationships of *MORF* genes in poplar. (**B**) Interchromosomal relationships of *MORF* genes in *A. thaliana*. Black lines indicated all synteny blocks in the Populus or *A. thaliana* genome and the red indicated synteny blocks where *MORF* duplicated gene pairs were. The chromosome number was indicated at the bottom of each chromosome and each cell on the outside of chromosomes represented 5 Mbp in *Populus* and 1.5 Mbp in *A. thaliana*. **Figure 5.** Interchromosomal relationships of *MORF* genes. (**A**) Interchromosomal relationships of *MORF* genes in poplar. (**B**) Interchromosomal relationships of *MORF* genes in *A. thaliana*. Black lines indicated all synteny blocks in the Populus or *A. thaliana* genome and the red indicated synteny blocks where *MORF* duplicated gene pairs were. The chromosome number was indicated at the bottom of each chromosome and each cell on the outside of chromosomes represented 5 Mbp in *Populus* and 1.5 Mbp in *A. thaliana*.

We further inferred the phylogenetic relationships of *PtrMORFs*, based on the diverse roles of DYW-PPR protein binding to MORFs on RNA editing sites in *A. thaliana* [9]. We constructed two syntenic maps of poplar with *A. thaliana* and rice. A total of four *PtrMORF* genes showed syntenic relationships with *MORFs* in *A. thaliana* (Figure 6A,C and Table S3). However, no synteny was in rice (Figure 6B). We performed phylogenetic analyses of these MORF proteins in Populus and *A. thaliana.*  The ratio of nonsynonymous/synonymous substitution rates of *PtrMORFs* and *MORFs* in *A. thaliana*  was determined to evaluate the selection pressure on amino acid substitutions (ω = dN/dS) and the role of Darwinian positive selection in driving gene divergence after duplication [35–37]. Generally, ω > 1 indicates positive selection, ω < 1 provides evidence for negative or purifying selection, and ω We further inferred the phylogenetic relationships of *PtrMORFs*, based on the diverse roles of DYW-PPR protein binding to MORFs on RNA editing sites in *A. thaliana* [9]. We constructed two syntenic maps of poplar with *A. thaliana* and rice. A total of four *PtrMORF* genes showed syntenic relationships with *MORFs* in *A. thaliana* (Figure 6A,C and Table S3). However, no synteny was in rice (Figure 6B). We performed phylogenetic analyses of these MORF proteins in Populus and *A. thaliana.*  The ratio of nonsynonymous/synonymous substitution rates of *PtrMORFs* and *MORFs* in *A. thaliana*  was determined to evaluate the selection pressure on amino acid substitutions (ω = dN/dS) and the role of Darwinian positive selection in driving gene divergence after duplication [35–37]. Generally, ω > 1 indicates positive selection, ω < 1 provides evidence for negative or purifying selection, and ω = 1 supports neutral evolution. We further inferred the phylogenetic relationships of *PtrMORFs*, based on the diverse roles of DYW-PPR protein binding to MORFs on RNA editing sites in *A. thaliana* [9]. We constructed two syntenic maps of poplar with *A. thaliana* and rice. A total of four *PtrMORF* genes showed syntenic relationships with *MORFs* in *A. thaliana* (Figure 6A,C and Table S3). However, no synteny was in rice (Figure 6B). We performed phylogenetic analyses of these MORF proteins in Populus and *A. thaliana*. The ratio of nonsynonymous/synonymous substitution rates of *PtrMORFs* and *MORFs* in *A. thaliana* was determined to evaluate the selection pressure on amino acid substitutions (ω = dN/dS) and the role of Darwinian positive selection in driving gene divergence after duplication [35–37]. Generally, ω > 1 indicates positive selection, ω < 1 provides evidence for negative or purifying selection, and ω = 1 supports neutral evolution.

**Figure 6.** Synteny analysis of MORF genes between poplar and two representative plant species. (**A**,**B**) Gray lines in the background indicate the collinear blocks within poplar and other plant genomes, while the red lines highlight the syntenic MORF gene pairs. (**C**) The phylogenetic relationship of 7 members of four pairs. (**D**) The Bayes Empirical Bayes (BEB) probabilities for sites in the positively selected class (ω > 1). The x-axis denotes position in the amino acid alignment. **Figure 6.** Synteny analysis of MORF genes between poplar and two representative plant species. (**A**,**B**) Gray lines in the background indicate the collinear blocks within poplar and other plant genomes, while the red lines highlight the syntenic MORF gene pairs. (**C**) The phylogenetic relationship of 7 members of four pairs. (**D**) The Bayes Empirical Bayes (BEB) probabilities for sites in the positively selected class (ω > 1). The x-axis denotes position in the amino acid alignment.

Using the maximum likelihood method and codon substitution models implemented in PAML, the selection pressure in the four groups of *MORF* genes was evaluated by likelihood ratio tests (LRTs). The estimated ω value for all MORF genes was 0.055 using a one-ratio model (M0). We then detected the positive selection acting on particular groups using a branch model in which each clade had its own ω value. The LRT statistic suggested that the ω values for groups II, III, and IV were significantly different from that of group I, and the ω estimates for all groups other than group I still suggested purifying selection (Table 1). Using the maximum likelihood method and codon substitution models implemented in PAML, the selection pressure in the four groups of *MORF* genes was evaluated by likelihood ratio tests (LRTs). The estimated ω value for all MORF genes was 0.055 using a one-ratio model (M0). We then detected the positive selection acting on particular groups using a branch model in which each clade had its own ω value. The LRT statistic suggested that the ω values for groups II, III, and IV were significantly different from that of group I, and the ω estimates for all groups other than group I still suggested purifying selection (Table 1).

Thus, we estimated the evolutionary forces acting on individual codon sites using site-specific likelihood models of codon substitution because positive selection was unlikely to affect all sites over prolonged time periods. Three pairs of models—M1 (neutral) and M2 (selection), M0 (one ratio) and M3 (discrete), and M7 (beta) and M8 (beta & ω)—formed three LRTs. As shown in Table 2, model M1 was not significantly worse than M2, although it suggested that 6.6% of sites were nearly neutral with ω = 1. Model M3 with K = 8 suggested that 0.7% of sites were under positive selection, and model M3 was significantly better than the one rate model. Model M8, in which an additional ω ratio was estimated from the data, was not significantly better than M7, indicating that no sites were under positive selection. Thus, we estimated the evolutionary forces acting on individual codon sites using site-specific likelihood models of codon substitution because positive selection was unlikely to affect all sites over prolonged time periods. Three pairs of models—M1 (neutral) and M2 (selection), M0 (one ratio) and M3 (discrete), and M7 (beta) and M8 (beta & ω)—formed three LRTs. As shown in Table 2, model M1 was not significantly worse than M2, although it suggested that 6.6% of sites were nearly neutral with ω = 1. Model M3 with K = 8 suggested that 0.7% of sites were under positive selection, and model M3 was significantly better than the one rate model. Model M8, in which an additional ω ratio was estimated from the data, was not significantly better than M7, indicating that no sites were under positive selection.

To better detect positive selection, the branch-site model was also applied to evaluate selection on all amino acids of MORF proteins in specific groups (Table 3). LRT showed that model A fitted the data significantly better than the site-specific model M1 (*p* < 0.01) in group II, implying positive selection on 29.7% of sites in group II. At the posterior probability (*p*) > 95%, four sites were likely to be under positive selection in group II. Referring to first sequence in group II, *MORF8*, these positively selected sites were 8T, 33R, 78D, and 83V (Figure 6D). However, no positive selection was found in the branch including *PtrMORFs*. To better detect positive selection, the branch-site model was also applied to evaluate selection on all amino acids of MORF proteins in specific groups (Table 3). LRT showed that model A fitted the data significantly better than the site-specific model M1 (*p* < 0.01) in group II, implying positive selection on 29.7% of sites in group II. At the posterior probability (*p*) > 95%, four sites were likely to be under positive selection in group II. Referring to first sequence in group II, *MORF8*, these positively selected sites were 8T, 33R, 78D, and 83V (Figure 6D). However, no positive selection was found in the branch including *PtrMORFs*.


**Table 1.**Parameter estimates and likelihood ratio tests for the branch model.


 **2.**Parameter estimates and likelihood ratio tests for the site model.

**Table**


**Table 3.** Parameter estimates and likelihood ratio tests for the branch-site model.

*Stress* 

#### *2.4. Expression of PtrMORF Genes and Six Genes from Chloroplasts and Mitochondria under Drought Stress* Although the role of MORF proteins in plastid development and RNA editing in *A. thaliana* and

*Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 10 of 19

Although the role of MORF proteins in plastid development and RNA editing in *A. thaliana* and rice have been studied, little is known about how *MORFs* respond to abiotic stimuli, such as drought stress, particularly in poplar [27,32]. rice have been studied, little is known about how *MORFs* respond to abiotic stimuli, such as drought stress, particularly in poplar [27,32]. Quantitative real-time reverse transcription-PCR (qRT–PCR) for long-term water deficiency

Quantitative real-time reverse transcription-PCR (qRT–PCR) for long-term water deficiency stress in black poplar (Populus × *euramericana cv. 'Neva'*) was further performed to evaluate the differential expression of each *PtrMORF* gene. Nine genes exhibited significantly different expression under limited water stress except *PtrMORF2.1*, and these could be preliminarily considered as drought-responsive genes. Among them, the two *PtrMORF* genes *PtrMORF1.3* and *PtrMORF8.1* were highly expressed followed by low expression. The expression levels of *PtrMORF1.1* and *PtrMORF1.2*, *PtrMORF2.2*, *PtrMORF3*, and *PtrMORF8.2* under drought conditions were significantly lower than those in control conditions. The mRNA accumulation of *MORF9* fluctuated obviously, exhibiting decreased expression after 3-day and 9-day drought, but increased approximately 1.5-fold after 6-day and 12-day drought. However, the pattern of the response to drought was not consistent. For example, the expression level of *PtrMORF1.1* had downregulated after 3-day drought stress. At the same time, *PtrMORF1.2* and *PtrMORF8.2* were significantly downregulated until 12 days of limited water stress. Furthermore, *PtrMORF* genes that belonged to the same group in the phylogenetic tree had different expression patterns under stress treatment. For instance, *PtrMORF1.1* and *PtrMORF1.3* were downregulated and upregulated in 3-day drought treatments, while *PtrMORF1.2*, also orthologous to group I members and *MORF1* in *A. thaliana*, was not affected by limiting water until 12 days post-treatment (Figure 7). stress in black poplar (Populus *× euramericana cv. 'Neva'*) was further performed to evaluate the differential expression of each *PtrMORF* gene. Nine genes exhibited significantly different expression under limited water stress except *PtrMORF2.1*, and these could be preliminarily considered as drought-responsive genes. Among them, the two *PtrMORF* genes *PtrMORF1.3* and *PtrMORF8.1* were highly expressed followed by low expression. The expression levels of *PtrMORF1.1* and *PtrMORF1.2*, *PtrMORF2.2*, *PtrMORF3*, and *PtrMORF8.2* under drought conditions were significantly lower than those in control conditions. The mRNA accumulation of *MORF9* fluctuated obviously, exhibiting decreased expression after 3-day and 9-day drought, but increased approximately 1.5-fold after 6-day and 12-day drought. However, the pattern of the response to drought was not consistent. For example, the expression level of *PtrMORF1.1* had downregulated after 3-day drought stress. At the same time, *PtrMORF1.2* and *PtrMORF8.2* were significantly downregulated until 12 days of limited water stress. Furthermore, *PtrMORF* genes that belonged to the same group in the phylogenetic tree had different expression patterns under stress treatment. For instance, *PtrMORF1.1* and *PtrMORF1.3* were downregulated and upregulated in 3-day drought treatments, while *PtrMORF1.2*, also orthologous to group I members and *MORF1* in *A. thaliana*, was not affected by limiting water until 12 days post-treatment (Figure 7).

**Figure 7.** Quantitative RT-PCR analysis of *PtrMORF* genes expression in response to drought stress. 0d, 3d, 6d, 9d, and 12d, drought for 0, 3 ,6, 9, and 12 days in a greenhouse environment, respectively. Data were normalized to UBQ gene. The sample without drought treatment was defined as 1 in the figure. The data were presented as the mean ± SE of three separate measurements. Asterisks denote significant differences: \* *p* < 0.05; \*\* *p* < 0.01. **Figure 7.** Quantitative RT-PCR analysis of *PtrMORF* genes expression in response to drought stress. 0d, 3d, 6d, 9d, and 12d, drought for 0, 3, 6, 9, and 12 days in a greenhouse environment, respectively. Data were normalized to UBQ gene. The sample without drought treatment was defined as 1 in the figure. The data were presented as the mean ± SE of three separate measurements. Asterisks denote significant differences: \* *p* < 0.05; \*\* *p* < 0.01.

In order to further explore the possible relationship between RNA editing and plant stress, we selected a total of six genes from chloroplasts and mitochondria to evaluate their expression, based on previous studies [27,38,39]. The *RPS14* gene represented important chloroplast proteins, and the In order to further explore the possible relationship between RNA editing and plant stress, we selected a total of six genes from chloroplasts and mitochondria to evaluate their expression, based on previous studies [27,38,39]. The *RPS14* gene represented important chloroplast proteins, and the fluctuant editing efficiency and gene expression of *RPS14* were associated with cytokinins

against stress [39,40]. Under salt stress, *PSBF* and *NDHB* likely showed an elevated editing percentage, linked to PSII repair and increase in *NDHB* gene translation [39]. In rice, *atp6*, *atp9*, and *ccmC* were related to mitochondrial electron transport chain and had been selected to test RNA editing during stress exposure [38]. In our study, these reported genes in poplar were obtained from National Center for Biotechnology information and showed different expression owing to drought. Six genes, except *PSBF*, were significantly upregulated or downregulated after drought stress (Figure 8). The three mitochondrial genes—*atp6*, *atp9*, and *ccmC*—were significantly downregulated after 3-day drought. The level of *RPS14* was highly expressed on 3-day drought followed by a low expression on 6-day drought, and the expression of poplar *NDHB* was significantly increased under 12-day drought. stress [39,40]. Under salt stress, *PSBF* and *NDHB* likely showed an elevated editing percentage, linked to PSII repair and increase in *NDHB* gene translation [39]. In rice, *atp6*, *atp9*, and *ccmC* were related to mitochondrial electron transport chain and had been selected to test RNA editing during stress exposure [38]. In our study, these reported genes in poplar were obtained from National Center for Biotechnology information and showed different expression owing to drought. Six genes, except *PSBF*, were significantly upregulated or downregulated after drought stress (Figure 8). The three mitochondrial genes—*atp6*, *atp9*, and *ccmC*—were significantly downregulated after 3-day drought. The level of *RPS14* was highly expressed on 3-day drought followed by a low expression on 6-day drought, and the expression of poplar *NDHB* was significantly increased under 12-day drought.

fluctuant editing efficiency and gene expression of *RPS14* were associated with cytokinins against

**Figure 8.** Expressions of three plastid (RPS14, PSBF, and NDHB) and three mitochondrial genes (atp6, atp9, and ccmC) in poplar on the 0, 3rd, 6th, 9th, and twelfth day after drought tolerance. Data were normalized to UBQ gene. The sample without drought treatment was defined as 1 in the figure. The data were presented as the mean ± SE of three separate measurements. Asterisks denote significant **Figure 8.** Expressions of three plastid (RPS14, PSBF, and NDHB) and three mitochondrial genes (atp6, atp9, and ccmC) in poplar on the 0, 3rd, 6th, 9th, and twelfth day after drought tolerance. Data were normalized to UBQ gene. The sample without drought treatment was defined as 1 in the figure. The data were presented as the mean ± SE of three separate measurements. Asterisks denote significant differences: \* *p* < 0.05; \*\* *p* < 0.01.
