Genome-Wide Identification and Evolutionary Analysis of AOMT Gene Family in Pomegranate (Punica granatum)

Gene duplication is the major resource with which to generate new genes, which provide raw material for novel functions evolution. Thus, to elucidate the gene family evolution after duplication events is of vital importance. Anthocyanin O-methyltransferases (AOMTs) have been recognized as being capable of anthocyanin methylation, which increases anthocyanin diversity and stability and improves the protection of plants from environmental stress. Meanwhile, no detailed identification or genome-wide analysis of the AOMT gene family members in pomegranate (Punica granatum) have been reported. Three published pomegranate genome sequences offer substantial resources with which to explore gene evolution based on the whole genome. Altogether, 58 identified OMTs from pomegranate and five other species were divided into the AOMT group and the OMT group, according to their phylogenetic tree and AOMTs derived from OMTs. AOMTs in the same subclade have a similar gene structure and protein conserved motifs. The PgAOMT family evolved and expanded primarily via whole-genome duplication (WGD) and tandem duplication. PgAOMTs expression pattern in peel and aril development by qRT-PCR verification indicated that PgAOMTs had tissue-specific patterns. The main fates of AOMTs were neoor non-functionalization after duplication events. High expression genes of PgOMT04 and PgOMT09 were speculated to contribute to “Taishanhong” pomegranate’s bright red peel color. Finally, we integrated the above analysis in order to infer the evolutionary scenario of AOMT family.


Introduction
Gene duplication and loss events constitute the main factor of gene family evolution [1]. Whole genome duplication (WGD) and small scale duplication events are the major sources of most duplicates [2]. After duplication events, duplicate copies are a predominant feature of individual genes and whole genomes, while a large fraction of duplications was lost [3]. Some of these duplicates' retention contributes to neo-functionalization, such as increasing metabolic activity, evolving new structure and novel adaptive traits [2,3]. After WGD, the strong biased retention of regulatory and developmental genes has ramifications for evolution in the longer term [4]. These regulatory and developmental genes control multiple reactions in land plants biochemical pathways. Thus, it is significant to elucidate the evolution of the regulatory and developmental gene families after WGD.
Anthocyanins, which are the largest group of flavonoids, are significant pigments to plants. By conferring to the fruit pulp, pericarp, flower and leaf, a red, orange or purple color, anthocyanin, contributes to their appealing appearance [5][6][7][8]. Through these bright colors, plants attract the coevolving birds or insects to spread seeds and pollens [9]. Anthocyanins have potential health benefits. Some research suggests that anthocyanins possess antioxidant, antitumor, anti-inflammatory and immunomodulation activity [10].
In order to verify the reliability of the AOMT candidate gene sequences, the results of RT-PCR and gene cloning were compared with the predicted genes. Fresh young leaves of pomegranate were selected as materials and stored at −80 • C after liquid nitrogen freezing.
In order to study the differential expression of AOMT genes in peel and aril at different stages during development, the fruits of four different development stages (July, August, September and October) were selected at random. The peel and aril were separated and stored at −80 • C after liquid nitrogen freezing.

Data Sources and Gene Identification of AOMTs
Pomegranate genome (ASM220158v1) was downloaded from NCBI (https://www. ncbi.nlm.nih.gov) [20]. A. thaliana AOMT protein sequences were downloaded from the A. thaliana database (http://www.arabidopsis.org). M. domestica, C. sinensis and V. vinifera AOMT protein sequences were downloaded from Phytozome and D. zibethinus was downloaded from NCBI. Detailed information of the species above was listed in the Appendix A.

Phylogenetic Analysis
To explore the phylogenetic relationship of the AOMT genes among six species, including 10 pomegranate OMT genes, 16 OMT D. zibethinus genes, 7 C. sinensis OMT genes, 6 A. thaliana OMT genes, 9 M. domestica OMT genes and 10 V. vinifera OMT genes, these identified AOMT gene sequences were utilized to perform multiple sequence alignment analysis with MAFFT (a multiple sequence alignment program) software [32]. Then, a maximum likelihood tree was established using PhyML [33] with the JTT model, SPR topology search, 1000 bootstrap replicates and aLRT statistics. The tree was visualized using Figtree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/).

Exon-Intron Structure
Gene Structure Displayer server (GSDS) website [34] (http://gsds.cbi.pku.edu.cn) was used to analyze the AOMT gene structure. Each plant species GFF annotation was downloaded corresponding to the genome. The required annotation contents were filtered out using perl script and were uploaded to the GSDS website for the AOMT exon-intron structure.

Motif Analysis
In order to further identified and analyzed the protein conserved motifs of the PgAOMT family members, The online tool Multiple Em for Motif Elicitation (MEME: http://meme.nbcr.net/meme/intro.html) [35] was employed with the following parameters: the motif width was 20~100 aa, and the number of motifs was 6.

RNA Extraction, cDNA Synthesis and RT-PCR Analysis
Total RNA of pomegranate young leaves was extracted using RNAprep Pure total RNA extraction kit. RNA was purified using DNase I digestion kit (Vazyme-innovation in enzyme technology, Nanjing, China). The obtained RNA was reverse transcribed into cDNA by HiScript ® II 1st Strand cDNA Synthesis Kit (Vazyme-innovation in enzyme technology, Nanjing, China). The pomegranate target genes were screened according to the verified A. thaliana, M. domestica, D. zibethinus, C. sinensis and V. vinifera AOMT homologous gene. Then, the primers were designed depending on the candidate genes with Primer5.0 for PCR amplification based on the predicted genes through ePCR and were listed in Appendix C. Finally, the GBclonart seamless cloning kit (GBI, Suzhou, China) was used to recombine gene as Zhang et al. [36] described and the cloned genes were sequenced.

Expression Analyses of PgAOMT Gene Family Members
The public transcriptomes pomegranate peel (SRR5279392, SRR5279393, SRR5279394) and aril (SRR5279386, SRR5279387, SRR5279388) during development were downloaded from NCBI. Quality controls for these transcriptomes were examined by the NGS QC Toolkit v2.3.3. Then, the gene expression was quantified by the kallisto v0.43.1. The data was calculated by reads per kilobase per million mapped reads (RPKM) as transcript abundance and the RPKM values were transformed in log2 fold change. The differently expressed PgAOMT genes in peel and aril across development periods were analyzed with DESeq2 software [37]. According to the threshold p-value adjust <0.01 and log2(fold-change) ≥2, differentially expressed genes (DEGs) were screened out. Up-regulated DEGs and downregulated DEGs were classified. The volcano map generation was exhibited by R software. To validate the result of the RNA-seq data, three samples for each tissue (peel and aril) were collected for the expression levels of PgAOMT genes.
RNA samples of aril and peel during development (July, August, September and October) were used for qRT-PCR experiments. The operation method of RNA extraction and cDNA reserve transcription is the same as 2.7 mentioned. qRT amplification was performed in triplicate using the Luna ® Universal qPCR Master Mix (NEB, Ipswich, MA, USA). The qRT-PCR primers of AOMTs listed in Appendix C. The expression levels of these genes were assessed using the method described by Liu et al. [38]. The correlative expression data was calculated according to the 2 −(∆∆CT) method [39]. Finally, the expression spectrum was draw out using R software.

Identification of the AOMT Genes in Pomegranate
To identify the AOMT genes of pomegranate, we used SelectHMM to (e-value ≤ 10) screen and CDD and SMART to identify all possible AOMT members which contain complete domains. Altogether 10 OMT candidate genes were identified from pomegranate genome.

Phylogenetic Analysis of PgAOMT Gene Family
A comprehensive and robust phylogenetic tree is the basis of studying phylogenetic genomics. To study the AOMT gene family evolutionary relationship, a total of 58 candidate OMT genes were identified from pomegranate (10 OMT genes), D. zibethinus (16 OMT genes), C. sinensis (7 OMT genes), A. thaliana (6 OMT genes), M. domestica (9 OMT genes) and V. vinifera (10 OMT genes) and constructed a maximum likelihood (ML) tree ( Figure 1) using PhyML [33]. This ML tree has robust bootstrap-supported values in general.

Phylogenetic Analysis of PgAOMT Gene Family
A comprehensive and robust phylogenetic tree is the basis of studying phylogenetic genomics. To study the AOMT gene family evolutionary relationship, a total of 58 candidate OMT genes were identified from pomegranate (10 OMT genes), D. zibethinus (16 OMT genes), C. sinensis (7 OMT genes), A. thaliana (6 OMT genes), M. domestica (9 OMT genes) and V. vinifera (10 OMT genes) and constructed a maximum likelihood (ML) tree ( Figure 1) using PhyML [33]. This ML tree has robust bootstrap-supported values in general.  Based on the phylogenetic tree branches and bootstrap-supported values, the OMT genes can be divided into Class I and Class II. In Class I, six species OMT genes evenly distributed, while in Class II, just D. zibethinus, M. domestica and V. vinifera were the main species. We defined Class I and Class II as AOMT group and OMT group, respectively ( Figure 1). Meanwhile, it was speculated that the AOMT group evolved from the OMT group.
It is apparent that each of the four species (A. thaliana, V. vinifera, M. domestica and D. zibethinus) constituted the gene cluster respectively in the basal part of the ML tree (Figure 1), indicating that OMT gene family has undergone species-specific expansion during evolution. In Class I, the branches were made up of GSVIVT01020603001, GSVIVT01015243001, GSVIVT01015246001, GSVIVT01015245001, GSVIVT01015244001, PgOMT06, PgOMT05, PgOMT08, PgOMT04, PgOMT03, PgOMT01 and PgOMT02, which represented a ratio of 5:7 of grape to pomegranate (Figure 1). It was indicated that pomegranate experienced one WGD event after paleohexaploidy shared with V. vinifera and considerable fraction duplicates were lost during the subsequent gene loss event [20,[41][42][43].

Gene Structure and Protein Conserved Motifs of PgAOMT Genes Family
To study the AOMT gene structure, we analyzed the distribution and amount of the exons and introns among six species OMT genes on the basis of phylogenetic similarities (Figure 2A,B). The results showed that the majority of OMT genes contained introns and the number of introns varied from 0 to 8 ( Figure 2B), only AT1G24733, MD02G1230800 and PgOMT08 had no intron. In the same subclade, OMTs seem to have similar intron/exons distribution, while there were some special cases such as GSVIVT01010469001, MD16G1118200 and PgOMT10, that were obviously different from the genes of the same subclade in length and number of intron (Figure 2A,B). In addition, the introns of early copy genes tended to be longer (Figure 2A,B). These supported the assumption that the length and number of introns facilitate gene family differentiation. Combined with previous phylogenetic analysis, it was speculated that gene duplication events with gene loss events influenced PgOMTs gene structure. in the first group were identical 100%. The alignment areas of the sequences of the remaining four groups were basically the same, although there were some differences. It was concluded that the five OMT genes are all AOMT genes and their quality was high, which can meet the requirements of the subsequent analysis. To research the distribution of AOMT proteins conserved motifs, six species protein sequences were analyzed using MEME. Most OMT proteins have six conserved motifs ( Figure 2C,D). However, some proteins such as AT1G24733 only had motif 1 and 2. Combined with the phylogenetic tree, it explained why AT1G24733 and AT1G24735 that formed an individual branch (Figure 2A,C). It was observed that four AOMTs including PgOMT01, PgOMT02, PgOMT03 and PgOMT04 belonging to the same subclade contained all 6 motifs. While motif distribution of three other AOMTs from the same subclade were different, as follows, PgOMT10 possesses motif 1, 2 and 6; PgOMT06 possesses motif 1 and 4; PgOMT05 contains motif 1, 2, 5 and 6 and PgOMT08 contains motif 3 and 6 ( Figure 2A,C). These illustrated that there were obvious differences among PgOMT proteins after the WGD event. The loss and retention of motifs may be related to the evolution of subclade.

RT-PCR Analysis of PgAOMT Candidate Gene Sequences
The RT-PCR technique was used to detect the reliability of five OMT genes. Five target sequences were obtained after reverse transcription. The sequencing reads of five OMT genes were aligned with predicted PgAOMTs downloaded from the database using MEGA-X software [44] for screening (Appendix D). Through alignment, it was found that although there were differences among several base pairs in partial regions, the structure of the five sequencing OMT genes was consistent with the predicted sequences on the whole.
The inconsistent regions (Appendix D: the red box areas) were randomly sampled. As can be seen from the enlarged partial drawing, the alignment between two sequences in the first group were identical 100%. The alignment areas of the sequences of the remaining four groups were basically the same, although there were some differences. It was concluded that the five OMT genes are all AOMT genes and their quality was high, which can meet the requirements of the subsequent analysis.

Expression Patterns of PgAOMT in Peel and Aril during Different Development
Some research has revealed that AOMT genes play roles in anthocyanin methylation that influence the plant coloration and improve self-protection from environmental stress [8,18]. To gain insight into the functions of the PgAOMT genes in the development of aril and peel, the expression patterns of PgOMT genes in aril and peel were analyzed using public RNA-Seq data from NCBI ( Figure 3A). The volcano plot analysis showed that PgOMT09 in aril was significantly down-regulated ( Figure 3A). The DEG analyses demonstrated that PgOMT09 played pivotal roles in aril development.
To further understand the physiological role of the PgAOMTs, we determined the PgAOMT gene family members expression in peel and aril during different developmental stages (July, August, September and October) ( Figure 3B). The histograms revealed that the transcript abundance of seven PgAOMT genes was different significantly in aril and peel, suggesting that there were differences in their functions of PgAOMT in aril and peel development. Notably, PgOMT04 exhibited high expression both in peel and aril development. PgOMT09 expression was down-regulated significantly at the early stage of aril development. PgOMT01, PgOMT02 and PgOMT03 belonging to the same subclade exhibited low or no expression. These phenomena indicated that during the genome duplication events, some genes (such as PgOMT09) evolved new functions and some genes (such as PgOMT01, PgOMT02 and PgOMT03) became nonfunctional pseudogenes due to segmental large-scale duplication events arose from tandem duplication [2]. Additionally, PgAOMT genes had aril/peel-specific expression patterns in pomegranate ( Figure 3B). For instance, PgOMT01 and PgOMT03 only expressed in aril. PgOMT05, PgOMT07 and PgOMT10 were more highly expressed in aril than in peel. exhibited low or no expression. These phenomena indicated that during the genome du-plication events, some genes (such as PgOMT09) evolved new functions and some genes (such as PgOMT01, PgOMT02 and PgOMT03) became nonfunctional pseudogenes due to segmental large-scale duplication events arose from tandem duplication [2]. Additionally, PgAOMT genes had aril/peel-specific expression patterns in pomegranate ( Figure 3B). For instance, PgOMT01 and PgOMT03 only expressed in aril. PgOMT05, PgOMT07 and PgOMT10 were more highly expressed in aril than in peel.

Discussion
AOMT enzymes are known to play roles in anthocyanin methylation. It can increase diversity and improve the stability of anthocyanin [8,18]. Few AOMTs were identified and characterized such as VvAOMT in grape, CkmOMT2 in cyclamen, AnthOMT in tomato, DIFe1, DIFe2a and DIFe2b in Petunia spp. and PsAOMT and PtAOMT in purple-flowered and red-flower Paeonia plants, respectively [8,17,18,45,46]. So far, few articles have described the evolution of the OMT gene family. Three published pomegranate genomes have been sequenced and these provide valuable resources to explore functional genes in genome-wide data [20,47,48]. In addition, pomegranate possesses diverse anthocyanins [49,50] and AOMTs expansion occurred after duplication events [20]. Here, we identified

Discussion
AOMT enzymes are known to play roles in anthocyanin methylation. It can increase diversity and improve the stability of anthocyanin [8,18]. Few AOMTs were identified and characterized such as VvAOMT in grape, CkmOMT2 in cyclamen, AnthOMT in tomato, DIFe1, DIFe2a and DIFe2b in Petunia spp. and PsAOMT and PtAOMT in purple-flowered and red-flower Paeonia plants, respectively [8,17,18,45,46]. So far, few articles have described the evolution of the OMT gene family. Three published pomegranate genomes have been sequenced and these provide valuable resources to explore functional genes in genome-wide data [20,47,48]. In addition, pomegranate possesses diverse anthocyanins [49,50] and AOMTs expansion occurred after duplication events [20]. Here, we identified and analyzed the expression patterns of AOMTs based on phylogenomic analyses to speculate on the evolution of PgAOMT gene family.
It has been reported that pomegranate shared the paleohexaploidy event with all eudicots and the paleotetraploidy event with E. grandis [20]. Gene duplication has been deemed as a predominant mechanism to increase functional diversification and the enhanced expression divergence of duplicated genes can substantially contribute to physiological diversification [51]. Hence, exploring the duplication events to illuminate gene functional diversification is significant. A total of 58 OMTs were identified from the selected six species. They were divided into Class I and Class II by constructing an ML tree that analyzed through bootstrap-supported values (Figure 1). Of the 10 candidate PgOMT genes, nine PgOMTs belong to Class I and only one belonged to Class II (Figure 1). We defined Class I as AOMT group and Class II as OMT group. It was speculated that AOMT group evolved from OMT group. In the subclade of Class I, there was a ratio of five grape to seven pomegranate, suggesting that there might be pomegranate genome duplication and loss events [20,52]. In contrast with the number of AOMTs in other species, PgAOMTs copy number was more than that of E. grandis (6) and C. sinensis (2) [53,54]. In our study, the number of AOMTs from pomegranate was more than that from other species (Figure 1), suggesting a lineage-specific expansion event occurred in pomegranate. High copy number of PgAOMTs may contribute to the diverse anthocyanins forming distinct colors of pomegranate [20].
By analyzing AOMT genes structure combined with the phylogenetic tree, it was obvious that the gene structure of AOMTs in Class I have some differences, but members in the same subclade seem to be similar (Figure 2A,B). Notably, in the clade of PgOMT06, PgOMT05 and PgOMT08, the three PgAOMTs have fewer introns than others ( Figure 2B). The absence of introns of the three PgAOMTs proved the presumption that there might be gene loss events in PgAOMT family due to neo-, sub-or non-functionalization [55,56]. Furthermore, we found that the motifs encoding PgAOMTs in earlier evolved clades were absent to varying degrees in Class I. For example, PgOMT10 only contained motif 1, 2 and 6; PgOMT06 only contained motif 1 and 4; PgOMT05 possessed motif 1, 2 and 5; PgOMT08 possessed motif 3 and 6, whereas the PgAOMTs that later evolved contained all 6 motifs (Figure 2A,C). In general, the AOMTs phylogenetic analysis is consistent with the gene structure.
Some gene families evolved and expanded through tandem duplication or segmental duplication along with high rates of birth and death. Gene family's expansion in these ways is of important significance to diversify their functions [57]. Vice versa, gene function may feedback on copy number and genome organization and thus give rise to the widely varying patterns of tandem or segmental duplication [57]. In brief, genes exist in the form of clusters after gene tandem duplication [58]. Deciphering gene clusters evolution is vital to offer new insight into gene family evolutionary scenario. By analyzing the gene mapping in pomegranate, it presented the existence of tandem duplication (Figure 1, PgOMT01, PgOMT02 and PgOMT03), which are identical to those confirmed by Yuan et al. [20]. Pomegranate genome large-scale duplication gave rise to AOMT tandem duplication [20]. Gene tandem duplication events drove the evolution of PgAOMT family to some extent.
Pseudogenes tend to have significantly lower expression compared with annotated genes [59]. They are considered as nonfunctional genes that have similar structure to functional genes in a domain family. Hence, numerous pseudogenes were misidentified as functional genes during the genome annotation process [59]. Three tandem duplicated genes of PgOMT01, PgOMT02 and PgOMT03 have low or no gene expression levels in peel or aril development ( Figure 3B), suggesting they were pseudogenes and genes may experience rapid birth and death during AOMT family evolution [60]. In addition, the expression of PgOMT04 and PgOMT10 demonstrated significantly higher than others both in peel and aril development ( Figure 3B). The high expression supported the assumption that the PgAOMTs evolved new functions after genome duplication [36].
Anthocyanin accumulation is known to be induced for resisting abiotic stresses [61]. AOMTs have been identified that are capable of anthocyanin methylation, particularly endows plant with an enhanced protection, have experienced genome duplication, presumably due to selection imposed by rapid changes in abiotic environments [18,20]. As we all know, methylation modification plays roles in increasing anthocyanin diversity and stability. In addition, the anthocyanin methylation level is significantly correlated with the hue towards red or purple of plant tissues [8,45,62]. Zhao et al. [63] proposed that AOMT could be responsible for pomegranate peel and aril color transition from white to red. It has been proposed that the methylation of B-ring hydroxyl groups generates the shift towards red [64]. In "Taishanhong" pomegranate, Cyanidin-3Glucoside and Delphinidin-3-glucoside were detected the major anthocyanins, but methylated anthocyanin contents have not been reported [65].
Subsequently, qRT-PCR ( Figure 3) results showed that AOMT genes were expressed in peel and aril of pomegranate fruit. PgOMT04 and PgOMT09 in peel and PgOMT04 in aril with higher expression throughout the development period, indicating that these O-methylations of anthocyanin and others may be closely related with fruit development and the ripening process. PgOMT09 was down-regulated continuously in aril during the early development process, suggesting that it played roles in early development. PgOMT03 and PgOMT10 in peel and aril and PgOMT10 in aril showed fluctuating expression pattern, suggesting that they may play roles in specific stages of fruit development. In contrast, PgOMT09 expressed both in peel and aril, with distinct expression pattern, indicating that there were different accumulation of O-methylation anthocyanin between them. Interestingly, PgOMT01 was specifically expressed in pomegranate peel, suggesting that its specific formation of O-methylated anthocyanin or others would occur in fruit peel.

Conclusions
In this study, we identified 10 OMTs from pomegranate genome and they were divided into two groups (AOMT and OMT group) by phylogenetic analysis. Pomegranate have undergone WGD, tandem duplication and gene loss events, presumably due to selection for improving their behavioral response to rapid changes of environment. The expression patterns of PgAOMTs indicated that some members function in anthocyanin methylation which enhances diversity and stability of anthocyanins and improves protection against stresses such as ultraviolet. This study also provides fundamental information on AOMTs for the coloration of pomegranate.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author and the public pomegranate transcriptomes presented in this study are available in insert article.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
A summary of species for phylogenetic analysis and their genomic source.

Appendix B
Simplified species tree of 6 species.