Comparative Analysis of Re-Annotated Genes Provides Insight into Evolutionary Divergence and Expressions of Aquaporin Family in Pepper

Aquaporins (AQPs) are known to have a vital role in water transport in all living organisms including agriculturally important crops, but a comprehensive genomic study of AQPs in pepper has not been implemented. Here, we updated previous gene annotations and generated a total of 259 AQP genes from five plants, including pepper. Phylogenetic and motif analyses revealed that a large proportion of pepper AQP genes belong to the specific subgroup of tonoplast intrinsic protein (TIP) subfamily, TIP4. Chromosomal localization and estimated duplication times illustrated that genes in TIP4 formed a tandem array on the short arm of chromosome 1, resulting from pepper-specific expansion after its divergence with Solanaceae species. Transcriptome analyses under various abiotic stress conditions revealed that transport-, photosystem-, and thylakoid-related genes were generally enriched in expression clusters containing AQP genes in pepper. These results provide valuable genomic resources and insight into the evolutionary mechanism that generate genomic diversity of the AQP gene family in pepper.


Introduction
Aquaporins (AQPs) are small transmembrane proteins (21-34 KDa) from a large family of major intrinsic proteins (MIPs) that facilitate the movement of water, glycerol, and small molecules such as carbon dioxide (CO 2 ) and hydrogen peroxide (H 2 O 2 ) through biological membranes [1,2]. In plants, the AQP gene family is responsible for the selective transport of water and regulates the response to abiotic stresses such as drought and saline conditions [3]. For example, AtPIP2;2 transports water in Arabidopsis [4], and SlTIP2;2 increases the transpiration of tomato under drought and salt stresses [5]. Aquaporins typically have six transmembrane α-helices (H1 to H6) that are connected by five spiral loops (LA-LE). Loops B and E have two conserved asparagine-proline-alanine (NPA) motifs. Groups of amino acid residues, two groups in helices 2 and 5 (H2 and H5) and two groups in loop E (LE1 and LE2), form the aromatic/Arg (ar/R) selectivity filter. The AQP gene family in plants is classified into five major subfamilies: plasma membrane intrinsic proteins (PIPs), tonoplast intrinsic proteins (TIPs), NOD26-like intrinsic proteins (NIPs), small basic intrinsic proteins (SIPs), and uncategorized X intrinsic proteins (XIPs) [6].
Hot pepper (Capsicum annuum) is an economically important crop belonging to the Solanaceae family with tomato and potato. It is not only a source of the popular spice capsaicin but is of nutritional value as a source of vitamins, minerals, and phenolic compounds [7][8][9]. According to the United Nations Food and Agricultural Organization (FAO, Rome, Italy), world production of hot pepper increased by about 28% over the last decade, and in 2019, the top 10 pepper-producing countries cultivated 32 million tons (http://www.fao.org/home/en/ (accessed on 2 December 2020)). Recently, pepper genome sequencing projects with transcriptome studies have enabled identification and characterization of important gene families related to agricultural traits [8,[10][11][12]; however, comprehensive understanding of AQP genes in pepper through integrative genomic and transcriptome studies has not been achieved.
In this study, we updated the annotation of AQP gene family including AQP-like as well as AQP genes in pepper, tomato, potato, rice, and Arabidopsis and performed a comparative analysis. Of the five species we studied, pepper had the most AQP genes, and 42% of these were in the TIP subfamily. Specifically, from phylogenetic analysis, we determined that most of these pepper genes (74%) were in the TIP4 subgroup. We also investigated the AQP domain motif composition of each subfamily of these five species. Motifs 7,19, and 5, located in downstream positions 11, 12, and 13, respectively, were specific to the TIP4 subgroup of pepper. We demonstrated, with chromosomal localization, that 52% of the pepper AQP genes in the TIP4 subgroup clustered in the short arm of chromosome 1. Most of these genes had been duplicated after the divergence time between pepper and Solanaceae species. Expression clustering and gene ontology (GO) enrichment analyses demonstrated that genes with putative functions related to chloroplast, thylakoid, and transport were enriched in clusters containing the pepper AQP genes, which may be evidence that these AQP genes work with other pepper genes in the plant response to abiotic stresses, especially in the transport system of thylakoid membranes.
We anticipate that these updated AQP gene models will provide novel genomic resources to horticulturalists and enrich understanding of recent evolutionary divergence in pepper AQP genes and the complex transcriptional network that operates in abiotic stresses.

Re-Annotation and Classification of AQP Genes in Pepper and Other Species
To reduce annotation bias from different annotation methods with resources [13], we performed re-annotation and generated a total of 259 AQP genes containing 29 (11.2%) newly annotated genes in addition to previously annotated genes (Table 1). In the updated AQP genes, we observed 73 pepper AQP genes that were 1.8-fold of AQP genes in rice and Arabidopsis genomes. In addition, numbers of AQP genes in other Solanaceae genomes were relatively higher than those of rice and Arabidopsis (Table 1). Table 1. Numbers of re-annotated AQP genes in five species. We classified AQP genes into five subfamilies by comparing specific amino acid sequences in ar/R selectivity filter positions among five genomes ( Table 2 and Table S2). A large number of the pepper AQP genes belonged to the TIP subfamily, approximately threefold more than in the other species. Specifically, 31 pepper AQP genes (42.5%) were in the TIP subfamily, suggesting that specific AQP genes have been expanded in the pepper genome (Table 1). We also observed many genes belong to the XIP subfamily in potato and tomato genomes, but none in Arabidopsis; this is consistent with a previous report [14,15]. These results demonstrated that the copy numbers of AQP genes in the five plant genomes have changed via proliferation of genes belonging to specific subfamilies in the Solanaceae, especially in pepper.  We investigated the motif composition of AQP domains to characterize the structural diversity of these domains, which are mostly conserved and essential for transport function of those genes ( Figure 1 and Figure S1, Table S3). Overall, motif compositions were similar in each subfamily among the five plant species. Specifically, the motifs at the N-and C-terminal regions were mostly conserved among subfamilies, whereas the central regions, which includes the subfamily-specific motifs with ar/R selectivity filter position, were relatively variable. Thus, it may be that the motifs in the central account for the structural diversity of AQP domains observed among subfamilies ( Figure 1 and Figure S1). When we examined the motifs of AQP domains in pepper, we observed the specific motifs (7,19,5) were enriched in the TIP subfamily compared to the other species ( Figure 1 and Figure S1). This suggests that those motifs contributed to the structural diversity of pepper AQP domains in the TIP subfamily.  Table S3. Numbers of motifs in pepper AQP domains are shown as a bar graph. Black arrows show locations of transmembrane helices and loops in motif positions of AQP domains. H1-H6, helix 1 to helix 6; LB, loop B; LE, loop E; PIP, plasma membrane intrinsic protein subfamily; TIP, tonoplast intrinsic protein subfamily; NIP, NOD26-like intrinsic protein subfamily; SIP, small basic intrinsic protein subfamily; XIP, uncategorized X intrinsic protein subfamily.

Phylogenetic Analysis of Plant AQP Genes
We next performed a phylogenetic analysis and determined the subgroups (of the subfamilies) to which the AQP genes belonged. The subfamilies PIP, TIP, and NIP were further separated into three, four, and three subgroups, respectively, whereas SIP and XIP remained as single groups ( Figure 2A). Notably, a large number of AQP pepper genes in the TIP subfamily were mostly in the TIP4 subgroup-23 genes, representing 74% of the pepper genes in the TIP subfamily ( Figure 2B). Examination of the motif composition of AQP domains in TIP subgroups revealed a set of specific motifs, 7, 9, and 15 on downstream regions, were only observed in pepper AQP domains of TIP4 ( Figure 2C). These results are consistent with proliferation of specific TIP genes in pepper genome via lineage-specific duplication, which led to expansion of pepper TIP genes ( Figure 2A).  Table S3. Numbers of motifs in pepper AQP domains are shown as a bar graph. Black arrows show locations of transmembrane helices and loops in motif positions of AQP domains. H1-H6, helix 1 to helix 6; LB, loop B; LE, loop E; PIP, plasma membrane intrinsic protein subfamily; TIP, tonoplast intrinsic protein subfamily; NIP, NOD26-like intrinsic protein subfamily; SIP, small basic intrinsic protein subfamily; XIP, uncategorized X intrinsic protein subfamily.

Phylogenetic Analysis of Plant AQP Genes
We next performed a phylogenetic analysis and determined the subgroups (of the subfamilies) to which the AQP genes belonged. The subfamilies PIP, TIP, and NIP were further separated into three, four, and three subgroups, respectively, whereas SIP and XIP remained as single groups (Figure 2A). Notably, a large number of AQP pepper genes in the TIP subfamily were mostly in the TIP4 subgroup-23 genes, representing 74% of the pepper genes in the TIP subfamily ( Figure 2B). Examination of the motif composition of AQP domains in TIP subgroups revealed a set of specific motifs, 7, 9, and 15 on downstream regions, were only observed in pepper AQP domains of TIP4 ( Figure 2C). These results are consistent with proliferation of specific TIP genes in pepper genome via lineage-specific duplication, which led to expansion of pepper TIP genes (Figure 2A).

Chromosomal Location and Gene Duplication Analyses of Plant AQP Genes
To identify the location of AQP genes, we mapped their physical positions on the chromosomes of each genome (Figures 3 and S2). Overall, AQP genes were distributed across all chromosomes, with a few exceptions. We mapped 64 AQP genes in pepper, which has 12 chromosomes, and 33% of them mapped to chromosome 1 ( Figure 3). More specifically, 12 of the pepper AQP genes in TIP4 consisted of a tandem array on the short arm of chromosome 1 (Figure 3). Considering that there were no TIP4 genes in the syntenic regions of the tomato or potato chromosomes, this suggests that these 12 pepper genes emerged recently as a result of pepper-specific gene duplication. We also noted that seven potato and five tomato AQP genes belonging to the XIP1 subgroup formed clusters on the long arm of chromosome 10, which was specific to the Solanum spp. (Figure S2).

Chromosomal Location and Gene Duplication Analyses of Plant AQP Genes
To identify the location of AQP genes, we mapped their physical positions on the chromosomes of each genome (Figure 3 and Figure S2). Overall, AQP genes were distributed across all chromosomes, with a few exceptions. We mapped 64 AQP genes in pepper, which has 12 chromosomes, and 33% of them mapped to chromosome 1 (Figure 3). More specifically, 12 of the pepper AQP genes in TIP4 consisted of a tandem array on the short arm of chromosome 1 (Figure 3). Considering that there were no TIP4 genes in the syntenic regions of the tomato or potato chromosomes, this suggests that these 12 pepper genes emerged recently as a result of pepper-specific gene duplication. We also noted that seven potato and five tomato AQP genes belonging to the XIP1 subgroup formed clusters on the long arm of chromosome 10, which was specific to the Solanum spp. (Figure S2).
We next explored the duplication history of AQP genes. Gene duplication time was estimated by calculating the synonymous substitution rates (Ks) between recently duplicated gene pairs (Figure 4). Our analysis showed that the pepper AQP genes, especially those in the TIP4 subgroup, were actively duplicated, with a maximum Ks value of approximately 0.15. This also suggests that most of the duplication of the pepper AQP genes occurred after the divergence between pepper and other Solanum species, which have a Ks value of 0.3 [8] (Figure 4B,C). Taken together, these results indicate that recent gene duplications mainly occurred in TIP4 subgroup was a major contributor for construction of diversified AQP gene repertories in the pepper genome through expansion of specific AQP genes. We next explored the duplication history of AQP genes. Gene duplication time was estimated by calculating the synonymous substitution rates (Ks) between recently duplicated gene pairs (Figure 4). Our analysis showed that the pepper AQP genes, especially those in the TIP4 subgroup, were actively duplicated, with a maximum Ks value of approximately 0.15. This also suggests that most of the duplication of the pepper AQP genes occurred after the divergence between pepper and other Solanum species, which have a Ks value of 0.3 [8] (Figure 4B,C). Taken together, these results indicate that recent gene duplications mainly occurred in TIP4 subgroup was a major contributor for construction of diversified AQP gene repertories in the pepper genome through expansion of specific AQP genes.

Expression Profiles and Clustering Analyses of Pepper AQP Genes from Plants Grown under Abiotic Stress
Plants adapt to stressful environments through various systemic physiological changes. We examined responses to abiotic stresses in pepper at the level of transcription and expression. We first established an expression profile for the AQP genes under cold, heat, salt, and mannitol stresses ( Figure S3). Overall, a few genes (14) belonging to the PIP subfamily were highly expressed in all stress conditions, whereas several genes (14) belonging to the TIP subfamily were not, with one exception: CaTIP4;10 was highly expressed in all stress treatments, suggesting that it may have a particular function in various abiotic stresses. To group pepper genes with parallel expression patterns and thus predict genes with similar functions, we then performed expression clustering analyses and identified 18,052 (50%), 17,438 (48.6%), 15,926 (44.3%), and 16,105 (45%) DEGs in plants subjected to cold, heat, salt, and mannitol stress, respectively. Among these DEGs, 36, 42, 35, and 35 were AQP genes in cold, heat, salt, and mannitol, respectively ( Figure  5A). The clusters containing the most AQP genes were the cold cluster 5 (36%, 13 genes),

Expression Profiles and Clustering Analyses of Pepper AQP Genes from Plants Grown under Abiotic Stress
Plants adapt to stressful environments through various systemic physiological changes. We examined responses to abiotic stresses in pepper at the level of transcription and expression. We first established an expression profile for the AQP genes under cold, heat, salt, and mannitol stresses ( Figure S3). Overall, a few genes (14) belonging to the PIP subfamily were highly expressed in all stress conditions, whereas several genes (14) belonging to the TIP subfamily were not, with one exception: CaTIP4;10 was highly expressed in all stress treatments, suggesting that it may have a particular function in various abiotic stresses. To group pepper genes with parallel expression patterns and thus predict genes with similar functions, we then performed expression clustering analyses and identified 18,052 (50%), 17,438 (48.6%), 15,926 (44.3%), and 16,105 (45%) DEGs in plants subjected to cold, heat, salt, and mannitol stress, respectively. Among these DEGs, 36, 42, 35, and 35 were AQP genes in cold, heat, salt, and mannitol, respectively ( Figure 5A). The clusters containing the most AQP genes were the cold cluster 5 (36%, 13 genes), the heat-cluster 4 (26%, 11 genes), the salt cluster 6 (29%, 10 genes), and the mannitol cluster 1 (26%, 9 genes) ( Figure 5B). Most of the differentially expressed AQPs were in the PIP subfamily, and all gene clusters from plants exposed to heat and mannitol treatments contained PIP genes. The plants subjected to heat stress had 15 genes in the TIP subfamily, more than plants exposed to the other stresses. Furthermore, functional enrichment analysis, performed to predict the putative gene functions, revealed that genes with chloroplast-, photosystem-, organelle-, and thylakoidrelated GO terms were generally enriched in the clusters including AQP genes, suggesting that those genes could be involved in response to various abiotic stresses along with AQP genes ( Figure 5C). Specifically, many plants use osmotic adjustment and ion transport through AQP [16] to adapt to stress, and the finding that anion and organic substance transport-related GO terms were enriched in the salt cluster 5 ( Figure 5C) supports this suggestion. Moreover, enrichment of GO terms related to the thylakoid membrane in cold Furthermore, functional enrichment analysis, performed to predict the putative gene functions, revealed that genes with chloroplast-, photosystem-, organelle-, and thylakoidrelated GO terms were generally enriched in the clusters including AQP genes, suggesting that those genes could be involved in response to various abiotic stresses along with AQP genes ( Figure 5C). Specifically, many plants use osmotic adjustment and ion transport through AQP [16] to adapt to stress, and the finding that anion and organic substance transport-related GO terms were enriched in the salt cluster 5 ( Figure 5C) supports this suggestion. Moreover, enrichment of GO terms related to the thylakoid membrane in cold cluster 5 were consistent with previous studies, which showed that the TIP genes in Arabidopsis control the volume of the thylakoid lumen during water or light stress, thereby creating an optimal photosynthetic environment [17]. Our results suggest that the pepper AQP genes also have an important role in molecular responses to changing conditions and act in concert with other pepper genes.

Discussion
AQP is an important gene family involved in water transport and other small molecules such as CO 2 and H 2 O 2 plants, but the complete understanding of the AQP gene family in pepper has not been accomplished. For in-depth studies, an accurate annotation process involving structural annotation to fully predict gene structures containing target domains within the genome sequence is essential to provide a stronger foundation for biological investigation of living organisms [18]. However, reports of missing functional genes [19] or incomplete gene models [13] continue to be registered, generating bias in downstream studies.
In this study, we re-annotated AQP genes in five plant genomes using the target gene family annotation approach [20] to update AQP gene information via detection of previously omitted AQP genes. In general, whole-gene annotation methods, which have been applied to various species including model plants, have pre-processing steps such as repeat-masking to efficiently reduce computing power and manual work [21,22]. However, these steps often generate omission of genes in annotations [13], which are genes that exist in the assembled genome but are obscured by repeat-masking or unidentified due to other problems. Target gene family annotation methods were developed to complement these issues, and we conducted re-annotation using TGFam-Finder, developed to identify all target genes in assembled genome via intensive annotation of genomic regions containing target domain(s) [20].
Prior studies have also reported improvement of annotations of model plants such as Arabidopsis and rice, as well as human genome annotation by identification of missing genes in previous annotations [13,20]. Of the updated genes, we newly identified AQP genes omitted in previous annotations including rice and Arabidopsis and especially in the pepper genome. The updated AQP genes could serve as genomic resources for further comparative and evolutionary analyses on the basis of their accurate copy numbers in the assembled genomes.
Comparison of the five updated AQP gene models revealed that pepper had the highest number of AQP genes. Many previous studies reported that different copy numbers of gene families between species are mainly derived from species-specific evolution and resulted in genetic diversification [23,24]. We also observed that most pepper AQP genes were contained in the TIP4 subgroup with pepper-specific motif sequences, suggesting that these genes were lineage-specifically expanded in pepper, thereby increasing the genetic diversity in the AQP repertories.
Chromosomal localization and gene duplication history analyses showed that these pepper-specific genes in the TIP4 subgroup formed a tandem array in the short arm of chromosome 1, and the genes mainly emerged after the speciation time between pepper and other Solanum spp. A previous study also reported a similar case in which certain genes in the pepper nucleotide-binding, leucine-rich repeat (NLRs) gene family formed a tandem array in a specific chromosome via massive duplication, indicating that it occurred specifically in pepper [25]. Although most of their roles are unknown, a previous study reported that one gene, Pvr4, emerged eight million years ago after speciation between pepper and tomato, gaining a specific function in disease resistance for various potyviruses and showing that a pepper-specific NLR gene had an important role [25]. Similarly, the AQP genes in TIP4, which have undergone species-specific duplication after the pepper speciation time, may also possess important characteristics of agricultural value.
We further inspected the forward genetic screen using the expression level of pepper genes to identify those that are potentially involved in responses to different abiotic stress conditions. The genes belonging to the PIP were highly expressed. This was consistent with many reports that PIP is known to function in a variety of crops; for example, in Acacia auriculiformis, the AaPIP1-2 gene eliminated ROS and played a positive role in the drop stress; in Populus, the overexpression of PIP1;1 induced more sensitivity of transgenic in osmotic conventions; and in rice, the downregulation of PIP2;1 inhibited root growth [26][27][28]. Expression clustering and GO enrichment analyses on the pepper DEGs under abiotic stresses revealed clusters including AQP were generally enriched in GO terms related to the chloroplast, photosystem, organelle, and thylakoid, suggesting that genes involved in these functions may work together with AQP in response to various stresses [29]. Similarly, we observed newly identified genes in various clusters with other pepper genes. For example, CANN_MIP.PGAv.1.6.scaffold855.1 and CANN_MIP.chr06.1 belonging to salt cluster1 and salt cluster3, respectively, could function together with other pepper genes in cellular processes such as cell communication and membrane-bounded organelles under a salt stress condition. Although the actual functions of TIP4 genes in pepper are difficult to elucidate, their potential functions were predicted on the basis of the expression profiling and clustering analyses. Most genes belonging to the TIP4 subgroup were not significantly expressed, which is consistent with other reports that gene expression of newly emerged genes is not actively detected [30]. However, we found one gene, CA.PGAv.1.6.scaffold1432.12, in mannitol cluster2 that had GO enrichment in ribosome-related tasks, suggesting that this gene may also play a role in plant physiological processes related to ribosomes.
The updated AQP genes will be valuable genomic resources for further functional and breeding research in agriculturally valuable crops. Moreover, our findings provide insights into the evolutionary mechanism and expressions of pepper AQP genes that contribute to increasing genetic diversity and unveiling potential functions, respectively.

Re-Annotation of AQP Genes in Five Plant Genomes
To re-annotate putative aquaporin (AQP) genes including AQP-like genes, we first collected genome resources of Arabidopsis thaliana [31], Oryza sativa [32], Solanum lycopersicum [33], Solanum tuberosum [34], and Capsicum annuum [8], including assembled genomes, annotated genes, and RNA-Seq data (Table S1). The re-annotation was performed with TGFam-Finder v1.20 considering parameters described previously [20]. We used tsv files generated from the functional annotation of proteins in the five genomes using InterproScan-5 [35] as "TSV_FOR_DOMAIN_IDENTIFICATION". PF00230 (MIP) from Pfam database was used as the "TARGET_DOMAIN_ID". Finally, we obtained re-annotated AQP genes as updated annotations for the five plant genomes, including previously annotated genes from downloaded annotations and genes that were newly annotated from the TGFam-Finder.
We assigned gene names for re-annotated AQP genes. For previously annotated genes, if names had been assigned in previous studies, as was the case for the AQP genes in potato and rice, we matched new gene name to the original names. However, genes in Arabidopsis and tomato that had been annotated in previous studies that used different annotation versions. To address this problem, we compared the previously annotated AQP genes in updated annotations to their original versions with BLASTP, and the gene names with the highest identity score were used. Pepper has not been previously studied; therefore, all 73 re-annotated genes were given new names. In addition, newly annotated genes of other species were also assigned new names (Table 1, Table 2, and Table S2).

Classification and Phylogenetic Analysis of AQP Genes
The re-annotated AQP genes of the five species were aligned using MAFFT v7.470 [36], and unnecessary gaps in the alignments were trimmed with trimAL v1.4 (-gappyout) [37]. We obtained information on the ar/R selectivity filter from each of the four positions in H2 (37th position in the aligned sequences), H5 (143rd), LE1 (152nd), and LE2 (158th), on the basis of the location of the conserved NPA motifs in the alignments ( Table 2 and Table S2). For Arabidopsis, rice, tomato, and potato, subfamilies were classified according to the predicted amino acids in each of four positions of the protein, as described previously [38][39][40][41]. We classified the subfamilies of pepper AQPs, however, by considering the known ar/R selectivity filter positions of other four species. We compared unclassified genes to the classified AQP genes using BLASTP with the best match score. From this, we identified members of all five subfamilies of AQP genes: PIP, TIP, NIP, SIP, and XIP in the five species (Table 1, Table 2, and  Table S2).

Motif Analysis of Subfamilies in Pepper Genome
To identify conserved sequences of AQP domains in five plant species, we used MEME suite v5.1.1 (-mod zoops-nmotifs 50-minw 10-maxw 50-objfun se-markov_order 0) (Table S3) [44]. Then, the motif position in AQP domains were predicted using MAST v5.1.1 [45]. We determined the general 14 position constructing AQP domains using the identified motif sequences by manual inspections, considering sequence homology. We predicted the number of transmembrane helices and position from AQP domain using TMHMM (Figure 1) [46]. Finally, we counted motif frequencies located in the general position and visualized them using ggplot2 [47] in R package.

Location of AQP Genes on Chromosomes in Five the Plant Genomes
Chromosomal distribution of the re-annotated AQP genes was visualized using MapChart v2.32 [48]. Gene names were marked on the map, and members of the same subgroup, determined by the phylogenetic analysis in this study, were identified by the same color.

Estimation of Duplication Time of AQP Genes
We estimated the divergence times of the AQP genes by first identifying duplicated gene pairs with DupGen_Finder pipeline [49]. Next, the gene pairs were aligned with PRANK (-f = fasta-codon) [50]. We calculated the synonymous substitution rates (Ks) of each gene pair with the KaKs_Calculator v2.0 (-m MYN) [51]. Values of Ks < 2 were visualized with ggplot 2 [47] in an R package.

Expressional Clustering and GO Enrichment Analyses in the Pepper Genome
To examine expression profiles of pepper AQP genes, we collected RNA-Seq data for leaf tissues under various abiotic stress (cold, heat, mannitol, and salt) at different times generated from a previous study [52]. We first performed quality trimming with CLC Assembly Cell (CLC Bio, Aarhus, Denmark) to clear sequences of poor quality. Next, we mapped the RNA-Seq to the genome assembly of pepper using HISAT 2 [53] (-dtax). The mapped reads were quantified and the "Fragment Per Kilobase of transcript per Million mapped reads" (FPKM) values for the pepper AQP genes were calculated with StringTie [54] (-e-B-G). The overall expression levels of the pepper AQP genes in different conditions were visualized with log2 (FPKM + 1) values using pheatmap v1.0.12 [55] software in an R package.
We performed expression clustering of the differentially expressed genes (DEGs) in the whole annotated pepper genome, including the newly annotated AQP genes. Expression clustering was calculated from log2 transformation of fold changes between the FPKM values measured in plants under stress conditions and unstressed control plants. The DEGs that had a p value ≤ 0.05 were clustered using the Mfuzz program [56] in an R package. The number of expression clusters was determined by k-means analysis. The GO terms for genes in each cluster were identified using Omicsbox (https://www.biobam.com/ omicsbox/ (accessed on 5 January 2021)), and enriched GO terms were determined by Fisher's exact test (false discovery rates corrected p-value ≤ 0.01).
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/plants10061039/s1, Figure S1: Numbers of motifs in AQP domains of five plant species. Figure  S2: Location of AQP genes on chromosomes for four plant species. Figure S3: Expression profiling of pepper AQP genes using RNA-Seq data under various abiotic stresses. Table S1: Genomic resources of five plant species used in this study. Table S2: Detailed information of individual AQP genes in four plant genomes. Table S3: Statistics of motifs identified from five plant AQP genes by MEME.