*2.8. Expression Profiles of RsAP2 Genes under Abiotic Stress*

Quantitative Real-Time Fluorescent Polymerase Chain Reaction (qRT-PCR) was conducted to clarify the expression level change in RsAP2 genes with cold, salt and drought stress treatments in closed buds with three biological replicates and two technical replicates. Twenty RsAP2 genes were selected as examples according to the conditions of being confirmed as AP2-like ethylene-responsive transcription factors with annotations of the Pfam, SwissProt and NCBI Non-Redundant (NR) protein databases and had relatively higher FPKM values at the closed bud stage, while their promoters contained cis-acting elements response to abiotic stress and abscisic acid. Some RsAP2 genes showed significant changes in expression levels in response to abiotic stress treatments.

Overall, most of the RsAP2 genes were sensitive to cold stress, as 16 RsAP2 genes had significantly upregulated relative expression with cold treatment for 4 and 24 h compared to the control treated at room temperature and shielded from light (Figure 9). The peak RsAP2-42, 66, 102, 106 and 119 expression levels occurred at 4 h of cold stress, which may indicate a fast response situation.

**Figure 9.** Relative expression level of 20 selected RsAP2 genes in response to cold (**A**), salt (**B**) and drought (**C**) treatments. Actin gene was taken as the internal reference gene. The relative expression levels are presented with mean and standard deviation of qRT-PCR. The significance test is presented with asterisks: \* *p* < 0.05, \*\* *p* < 0.01, Student's *t*-test.

Although salt treatment may also cause drought, several RsAP2 genes, such as RsAP2- 30 and RsAP2-48, had obvious differences in expression patterns. RsAP2-30 expression had barely any volatility with the drought treatment, but the salt treatment resulted in this gene being significantly downregulated at 4 h but significantly upregulated at 24 h. As for RsAP2-30, the expression was mostly stable with the salt treatment, but the drought treatment resulted in extreme upregulation of the expression at 4 and 24 h. However, the expression patterns were rather similar in RsAP2-2, 15, 19, 40, 66, 106 and 116 with both the salt and drought treatments.

#### **3. Discussion**

With a large number of family members and diversified functions including abiotic stress resistance, hormone signal transduction and plant metabolite regulation, the AP2/ERF transcription family is important in many plant species. Many species specific AP2/ERF transcripts have been identified with the study of whole-genome sequencing data of many plants such as *Arabidopsis* [32], rice [33] and poplar [34]. Most studies have concentrated on model and crop plants, while only a few studies have detected transcription factor families in ornamental plants. Therefore, we performed a comprehensive and overall study of the AP2/ERF transcript family in *Rhododendron*, which has not been done before.

We identified 120 RsAP2 genes and carried out comprehensive studies on their phylogeny, sequence structure, chromosomal location, collinearity and expression patterns in flower development stages and abiotic stresses. This information can provide a theoretical basis and data support for further research on the functional pathways in stress resistance and flowering physiology of *Rhododendron* and other ornamental plants.

Previous research has reported that the number of AP2/ERF genes was not directly related to genome size, as the numbers of AP2/ERF genes in *Arabidopsis*, rice, poplar and rose were 147, 164, 200 and 135, respectively, while their genome sizes were 125, 430, 485 and 600 Mb, respectively. With 120 AP2/ERF genes identified and a 609 Mb genome size, the findings for *Rhododendron* could support this conclusion.

Phylogenetic trees can be used to study visually the relationships of sequences. In our research, the sequences were not clustered by species firstly but, rather, according to the classification of structure. The AP2 subfamily genes were grouped into an independent cluster, while the DREB and ERF subfamilies had a staggered distribution and formed different subgroups. This may be due to the tandem double AP2 domains, a relatively obvious specific structure.

The conserved motifs are conserved amino acid sequences with various biological functions involved in biological processes such as protein cointeraction, transcriptional activation and nuclear localization [32]. The tandem motif 11-1-8 and motif 3-2 occurred in all of the ERF and DREB subfamily members, while in AP2 subfamily members the tandem motif 6-4-5 was a common structure, suggesting that sequence structures have a close relationship with certain functions.

Some characteristics and commonalities can be found through the analysis of gene structures. There were 79 RsAP2 genes with no introns, which accounted for 65.83% of all the family members, similar to the situation in rose [7]. Interestingly, the RAV subfamily members all had no introns, while the AP2 subfamily members contained more than four introns without exception.

According to many studies, the evolution of genomes and genetic systems cannot occur without gene duplications. Gene families could therefore be the product of segmental duplication and tandem duplication [35]; in these ways, plants can achieve rapid adaption to environmental change [36]. In the present study, the duplication event of RsAP2 genes was analyzed to complete the understanding of the expansion pattern. The distribution of tandem duplications was uneven, with 30 RsAP2 genes clustered into 12 tandem duplication event regions on *Rhododendron* chromosomes 1, 5, 7, 8, 9 and 10, while 50 segmental duplication events involving 71 RsAP2 genes were identified. This indicated that tandem duplication and segmental duplication both contributed greatly to the evolution of *Rhodo-* *dendron* AP2/ERF genes. Many of the RsAP2 genes associated with tandem duplication events tended to have similar expression patterns in different developing stages of *Rhododendron* flowers, as shown in Figure 7, which is evidence of the close relationship between structure and function.

Cis-acting elements not only determine the tissue-specific expression of genes, but also have close relationship with stress resistance [37]. The cis-acting elements were mainly distributed in the 2 kb upstream sequence and could be analyzed to clarify the mechanisms and potential functions of RsAP2 transcriptional regulation. Among the elements, ABREs had the largest variety of sequence types according to data for four species (Table S2). The ABREs in *Rhododendron* were widely distributed upstream of 108 RsAP2 genes, and 78 of them appeared more than once. The high frequency and wide range of ABRE elements suggest that the function involved in abscisic acid responsiveness to resist abiotic stresses including cold, drought and salt [38], played an important role and was guaranteed in quantity. Another important cis-acting element is mainly involved in low temperature responsiveness [39], which occurred in the promoter of 68 RsAP2 genes. Three kinds of elements with MYB binding sites for MYB TFs, MBS, MBSI and MRE, were detected for their functions of drought-inducibility, flavonoid biosynthetic genes regulation and light responsiveness, respectively. This may imply the relation and interaction of AP2/ERF and MYB transcription factors and the complex regulatory network.

The expression patterns of genes were analyzed to achieve deeper understanding of the AP2/ERF genes in terms of overall trends and individual differences. Transcriptome data were the main data source to reveal the expression level of genes in different states of plant tissue. There is a period of dormancy with special biological significance for the *Rhododendron* flower bud before anthesis. According to the RNA-seq data for *R. simsii*, the RsAP2 genes presented different expression patterns in the five developmental stages of flowers. At the closed buds stage (S1), several RsAP2 genes had peak expression levels, while others had the lowest expression levels. Additionally, combined with the gene duplication events mentioned above, the clustered tandem duplicated genes exhibited high similarities in expression patterns among all five stages. This phenomenon may indicate the important roles of duplicated genes in developmental physiology and stress response because their function may be enhanced through an increase in quantity.

Quantitative RT-PCR can intuitively illustrate the RsAP2 genes in *Rhododendron* flower buds under different abiotic stress treatments. Twenty RsAP2 genes were selected as examples according to conditions, including phylogenetic analysis, annotations, FPKM values, cis-acting elements and protein interaction networks. With cold hardiness being an important breeding goal for researchers around the world, we focused on the gene response to cold stress, which significantly affected most of the selected RsAP2 genes. The accumulation of expression levels was positively correlated with the increase in time treated with cold stress for the overall trend. The functions of RsAP2-30/54/74/76 were supposed to be similar to those of the AT4G27950 gene, with a close phylogenetic relationship. AT4G27950 was regarded as CRF4 and showed transcriptional induction after exposure to cold [40]. Several special cases, such as RsAP2-42/106/119, had their expression peak with 4 h of cold treatment, which may indicate the characteristics of a rapid response. Although the application of higher salt concentrations may cause cellular osmotic water loss resulting in a drought-like effect, some AP2/ERF genes show specific patterns in response to salt treatment compared with mannitol-induced drought treatment, such as RsAP2-70 and RsAP2-48. RsAP2-70 had a close structural relationship with AT1G06160, which was named AtERF15. AtERF15 was reported to be the binding factor of Coupling Element 1 and was a positive regulator of the ABA response [41]. It can be assumed from this that RsAP2-70 was induced with salt stress and regulated by ABA as a downstream response factor. Meanwhile, RsAP2-48 exhibited no significant change in expression level after salt treatment. The function of RsAP2-48 may be similar to that of AT3G20840, the PLETHORA 1 gene in *Arabidopsis*, induced by auxin to control the alteration of stomatal apertures under drought stress [42]. Abiotic stresses, plant growth regulators and TFs

can form regulatory networks, and the AP2 gene family has extensive involvement in the process of environmental adaptation. Some of the RsAP2 genes appear to be redundant under abiotic stress and show similar expression patterns; they can be multiple types of insurance in terms of mechanisms or can be functionally cumulative in terms of quantity.

The above analysis was performed to systematically mine the genome-wide AP2 gene family of *Rhododendron*. Precise studies of regulatory systems and molecular mechanisms can provide a theoretical basis for future genetic improvement.

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

## *4.1. Identification of AP2 Family Members in the Rhododendron Genome*

The whole *Rhododendron* protein sequence was generated from the *Rhododendron* genome sequencing data. (The data have been submitted to NCBI genebank under the BioProject number of PRJNA902133). HMM was used as the main method to identify *Rhododendron* AP2 candidates. We downloaded the model file (PF00847) from the Protein Family (Pfam) database (http://pfam.xfam.org/, accessed on 23 June 2022). The first domain search was conducted using the model file with the standard of E-value < 1 × <sup>10</sup>−<sup>10</sup> to achieve a species-specific model in the *Rhododendron* genome protein sequence. Putative protein sequences were achieved with the second search using the species-specific model with the standard of E-value < 1 × <sup>10</sup>−10. The putative sequences were submitted to the SMART (http://smart.embl-heidelberg.de/, accessed on 28 June 2022), NCBI CDD (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi, accessed on 28 June 2022) and Pfam databases to verify the AP2 domain. The resulting sequences with AP2 domains confirmed by all three databases were clarified by removing redundant sequences, and the other sequences were categorized as RsAP2 genes. The RsAP2 genes were sorted according to their chromosome number and their location on the chromosome, and they were named in this order at the same time. The subcellular localization of RsAP2 proteins was predicted with WoLF PSORT (https://wolfpsort.hgc.jp/, accessed on 12 October 2022). The RsAP2 sequences were submitted to the Expert Protein Analysis System (ExPASy, http://web.expasy.org/protparam/, accessed on 28 June 2022) to analyze their molecular weights and theoretical isoelectric points (pI).
