Genome-Wide Identification and Function of Aquaporin Genes During Dormancy and Sprouting Periods of Kernel-Using Apricot (Prunus armeniaca L.)

Aquaporins (AQPs) are essential channel proteins that play a major role in plant growth and development, regulate plant water homeostasis, and transport uncharged solutes across biological membranes. In this study, 33 AQP genes were systematically identified from the kernel-using apricot (Prunus armeniaca L.) genome and divided into five subfamilies based on phylogenetic analyses. A total of 14 collinear blocks containing AQP genes between P. armeniaca and Arabidopsis thaliana were identified by synteny analysis, and 30 collinear blocks were identified between P. armeniaca and P. persica. Gene structure and conserved functional motif analyses indicated that the PaAQPs exhibit a conserved exon-intron pattern and that conserved motifs are present within members of each subfamily. Physiological mechanism prediction based on the aromatic/arginine selectivity filter, Froger’s positions, and three-dimensional (3D) protein model construction revealed marked differences in substrate specificity between the members of the five subfamilies of PaAQPs. Promoter analysis of the PaAQP genes for conserved regulatory elements suggested a greater abundance of cis-elements involved in light, hormone, and stress responses, which may reflect the differences in expression patterns of PaAQPs and their various functions associated with plant development and abiotic stress responses. Gene expression patterns of PaAQPs showed that PaPIP1-3, PaPIP2-1, and PaTIP1-1 were highly expressed in flower buds during the dormancy and sprouting stages of P. armeniaca. A PaAQP coexpression network showed that PaAQPs were coexpressed with 14 cold resistance genes and with 16 cold stress-associated genes. The expression pattern of 70% of the PaAQPs coexpressed with cold stress resistance genes was consistent with the four periods [Physiological dormancy (PD), ecological dormancy (ED), sprouting period (SP), and germination stage (GS)] of flower buds of P. armeniaca. Detection of the transient expression of GFP-tagged PaPIP1-1, PaPIP2-3, PaSIP1-3, PaXIP1-2, PaNIP6-1, and PaTIP1-1 revealed that the fusion proteins localized to the plasma membrane. Predictions of an A. thaliana ortholog-based protein–protein interaction network indicated that PaAQP proteins had complex relationships with the cold tolerance pathway, PaNIP6-1 could interact with WRKY6, PaTIP1-1 could interact with TSPO, and PaPIP2-1 could interact with ATHATPLC1G. Interestingly, overexpression of PaPIP1-3 and PaTIP1-1 increased the cold tolerance of and protein accumulation in yeast. Compared with wild-type plants, PaPIP1-3- and PaTIP1-1-overexpressing (OE) Arabidopsis plants exhibited greater tolerance to cold stress, as evidenced by better growth and greater antioxidative enzyme activities. Overall, our study provides insights into the interaction networks, expression patterns, and functional analysis of PaAQP genes in P. armeniaca L. and contributes to the further functional characterization of PaAQPs.

Aquaporins (AQPs) are essential channel proteins that play a major role in plant growth and development, regulate plant water homeostasis, and transport uncharged solutes across biological membranes. In this study, 33 AQP genes were systematically identified from the kernel-using apricot (Prunus armeniaca L.) genome and divided into five subfamilies based on phylogenetic analyses. A total of 14 collinear blocks containing AQP genes between P. armeniaca and Arabidopsis thaliana were identified by synteny analysis, and 30 collinear blocks were identified between P. armeniaca and P. persica. Gene structure and conserved functional motif analyses indicated that the PaAQPs exhibit a conserved exon-intron pattern and that conserved motifs are present within members of each subfamily. Physiological mechanism prediction based on the aromatic/arginine selectivity filter, Froger's positions, and three-dimensional (3D) protein model construction revealed marked differences in substrate specificity between the members of the five subfamilies of PaAQPs. Promoter analysis of the PaAQP genes for conserved regulatory elements suggested a greater abundance of cis-elements involved in light, hormone, and stress responses, which may reflect the differences in expression patterns of PaAQPs and their various functions associated with plant development and abiotic stress responses. Gene expression patterns of PaAQPs showed that PaPIP1-3, PaPIP2-1, and PaTIP1-1 were highly expressed in flower buds during the dormancy and sprouting stages of P. armeniaca. A PaAQP coexpression network showed that PaAQPs were coexpressed with 14 cold resistance genes and with 16 cold stressassociated genes. The expression pattern of 70% of the PaAQPs coexpressed with cold stress resistance genes was consistent with the four periods [Physiological dormancy (PD), ecological dormancy (ED), sprouting period (SP), and germination stage (GS)] of flower buds of P. armeniaca. Detection of the transient expression of GFP-tagged PaPIP1-1, PaPIP2-3, PaSIP1-3, PaXIP1-2, PaNIP6-1, and PaTIP1-1 revealed that INTRODUCTION Aquaporins (AQPs) are involved in water transport and facilitate the passage of other small solutes, such as CO 2 , boric acid, H 2 O 2 , glycerol, urea, and silicic acid, through cell membranes (Maurel et al., 2008(Maurel et al., , 2009. In plants, AQPs play a central role in maintaining the hydraulic conductivity balance and water homeostasis and are involved in the tolerance of abiotic stresses such as cold, drought, or salt stresses. In addition, the roles of AQPs in various growth processes, such as fruit ripening, petal and stomatal regulation, male fertility, seed germination, petal movement, and guard cell closure, have been documented (Heinen and Chaumont, 2009;Mitani-Ueno et al., 2011).
Kernel-using apricot (P. armeniaca L., syn. Armeniaca vulgaris Lam., also considered by some researchers to be P. armeniaca L. × P. sibirica L.) is an economically and ecologically important forest tree species in the three northern regions of China (northwestern China, northern China, and northeastern China). Because of its broad adaptability, strong resistance to drought, and easy management, it is one of the pioneer tree species used for ecological construction in mountainous areas and sandy lands, especially in western China. It is one of the six major nut-producing species in the world and is an important species for the production of woody oil and protein beverages. It can also contribute to the economies of poverty-stricken mountainous areas in China. The first high-quality apricot plant (Chuanzhihong) genome has been sequenced and released, laying a foundation for research on the molecular mechanisms underlying important agronomic characteristics of apricot (Jiang et al., 2019). Owing to its short dormancy period and early flowering period, kernel-using apricot (P. armeniaca L.) easily suffers from late frost damage, resulting in yield reduction or even no production. An effective strategy to solve the problem of a low and unstable yield of kernel-using apricot is to mine key genes related to flower bud dormancy and cold resistance, determine gene expression patterns during dormancy and sprouting periods (SPs), elucidate possible functional mechanisms and determine how the flowering period can be artificially manipulated.
Physiological dormancy (PD) (which is generally called endodormancy) is defined as the stage in which buds cannot burst even under optimal environmental conditions, while ecological dormancy (ED) (which is generally called ecodormancy) is the stage in which buds cannot burst because of lacking environmental factors (especially the temperature requirement accumulation for flower buds in ED periods is not enough). Endodormant buds, which cannot burst under preferred conditions due to internal physiological factors (such as embryos that are not mature physiologically, have a hard seed coat with poor permeability, or contain substances inhibiting germination) (Mitrakos and Diamantoglou, 2006;Gao et al., 2020), transition into the ecodormant (ecologically dormant) stage and acquire bud break capacity after exposure to a certain period of low temperature. When external environmental conditions (the temperature requirement accumulation for flower buds in ecologically dormant periods is enough) (Mazzitelli et al., 2007;Rohde and Bhalerao, 2007) are satisfied, the flower buds enter the SP and germination stage (GS), showing growth signs and sprouting, after which they can fully develop.
Cold stress can cause chilling injury in late spring, which obviously limits plant development and influences fruit yield in economic forests. Previous studies have shown that AQP family members play vital roles in different dormancy periods and lead to freezing tolerance and cold acclimation. For example, in A. thaliana, the expression of AtPIP2-5 and AtPIP2-6 is significantly upregulated by cold stress in roots and other parts of the plant (Alexandersson et al., 2010). The expression of rice OsTIP1-1 is downregulated under cold stress conditions (Sakurai et al., 2005) but upregulated in response to water and salinity stress (Liu et al., 1994). Transgenic studies on AQPs have been conducted to examine their responses to abiotic stresses such as cold, drought, and salinity. A. thaliana plants overexpressing AtPIP1;4 or AtPIP2-5 exhibit greater tolerance to cold stress (Afzal et al., 2016). Compared with wild-type tobacco plants, transgenic tobacco plants overexpressing the wheat TaAQP7 (PIP2) gene display improved cold tolerance (Zhou et al., 2012). Overexpression of MusaPIP1-2 in transgenic banana plants enhances tolerance to cold and drought conditions (Shekhawat and Ganapathi, 2013). Although AQPs have been extensively studied in both dicots and monocots through genome-wide and functional analyses (Nguyen et al., 2013), our knowledge of AQPs in P. armeniaca is still limited, and the expression profile of AQP members in flower buds in the different dormancy and sprouting stages is largely unknown in this species.
In this study, the genome sequence of P. armeniaca (Longwangmao) and the sequences of known AQPs from other species (such as A. thaliana and Prunus persica) were used to identify the genes encoding AQP members. The phylogenetic relationships, chromosome distribution, conserved residues and elements, protein-protein interactions, and expression patterns of AQPs were investigated to identify their functions in cold stress in kernel-using apricot. These results provide valuable information for further studies on the molecular mechanism of AQPs, which will facilitate the selection of candidate genes for cold tolerance improvement through genetic engineering in kernel-using apricot.

Identification of Aquaporin Family Members in P. armeniaca
The predicted peptide sequences of AQP members were acquired from the genome database (unpublished) of kernelusing apricot (P. armeniaca L.), also known as Longwangmao (Chinese Pinyin name), to construct a local protein database. BLASTP searches were performed using AQP protein sequences of A. thaliana and poplar (Table 1) as queries in The Arabidopsis Information Resource (TAIR) 1 and the Phytozome database (release version 12.0 2 ) as previously described (Gupta and Sankararamakrishnan, 2009;Zhang et al., 2013), with an E-value of 1e −10 and a minimum amino acid identity of 50%. The hidden Markov model (HMM) data of the MIP domain (PF00230) were downloaded from the Sanger database 3 . The PF00230 sequence information was then used to query the P. armeniaca protein database via HMMER 3.0 software 4 with default parameters. Thereafter, the predicted AQP protein sequences were submitted to the National Center for Biotechnology Information (NCBI) Conserved Domain Database (CDD) 5 (Marchler-Bauer et al., 2002 and SMART (Simple Modular Architecture Research Tool 6 ) (Letunic et al., 2012) to confirm the presence and completeness of the MIP domain, with an E-value threshold of 1e −2 .

Sequence Alignments and Phylogenetic and Synthetic Analyses
The PaAQP sequences were aligned with those of AtAQPs from A. thaliana 7 (Johanson et al., 2001) and those of PpeAQPs from P. persica (Deshmukh et al., 2015) via the ClustalW program 8 (Thompson et al., 1994). A phylogenetic tree was constructed by Molecular Evolutionary Genetics Analysis (MEGA) software 7.0 (Kumar et al., 2016;Qu et al., 2018) with the maximum likelihood method and 1,000 bootstrap resamplings. PaAQPs were named and confirmed by sequence homology and phylogenetic analyses. The syntenic relationships between the P. armeniaca, A. thaliana, and P. persica genomes were constructed by MCScan (Tang et al., 2008a(Tang et al., ,b, 2009 following the protocol for the Plant Genome Duplication Database (Lee et al., 2012). To identify the syntenic regions among P. armeniaca, initially, potential homologous gene pairs were detected using the BLASTP algorithm (E < 1e −5 , top five matches) between every two genomes. Second, the BLASTP results and the gene location information were input into MCScanX software (Wang et al., 2012b). P. armeniaca whole genome/segmental and tandem duplications were subsequently identified. All syntenic relationship and gene location data were presented using TBtools software .

Transient Expression of PaAQPs in A. thaliana Protoplasts
The GFP DNA fragment was amplified by PCR and subcloned into the pBI221 vector (Clontech, Mountain View, CA, United States). The coding regions of PaPIP1-1, PaPIP2-3, and PaSIP1-3 cDNA were amplified and fused in frame to the 5 end of the GFP DNA sequence. The six PaAQP-GFP constructs and the GFP-only control construct were introduced into A. thaliana protoplasts via a polyethylene glycol-mediated method described previously with some modifications (Sheen, 2001). For plasma membrane localization, the membrane-specific dye FM4-64 was used to stain the transfected protoplasts (Nelson et al., 2007). The fluorescence of FM4-64 and GFP in the cells was analyzed with a 551-nm helium-neon laser and a 488-nm argon laser, respectively, via confocal scanning microscopy (TCS SP8, Leica, Germany; GFP excitation at 488 nm and emission at 505-525 nm; FM4-64 excitation at 551 nm and emission at 498-588 nm).
Exon-Intron Structure, 3D Structure, and Conserved Motif Distribution The exon-intron organization was assessed via the Gene Structure Display Server 13 based on annotated coding sequences and genome sequences. Conserved domains and motifs were analyzed with MEME 14 . The three-dimensional (3D) structures of the PaAQPs were generated from protein sequences submitted to the Phyre2 server 15 . The transmembrane helices and topology of the PaAQPs were analyzed via the MEMSAT-SVM prediction method on the Phyre2 server.

In silico Analysis of Promoter Sequences and Protein-Protein Interaction Analysis
The promoter sequences (2.0 kb upstream region from the transcription start site) were extracted from the genomic DNA sequences of the PaAQP genes. Putative cis-acting regulatory elements of the promoter sequences were identified with the program PlantCARE online 16 . The prediction of interacting networks of proteins was generated from the STRING database 17 .

cDNA-Library Preparation and Illumina Sequencing for Transcriptome Analysis
Total RNA was extracted from the ED, SP, and GS of the flower buds of two Youyi plants using a Plant RNA Extraction Kit (Autolab, China) following the manufacturer's protocol. The concentration and quality of each RNA sample were determined using a NanoDrop 2000 TM microvolume spectrophotometer (Thermo Scientific, Waltham, MA, United States) and gel electrophoresis. Poly(A) mRNA was isolated from 10 mg of total RNA using magnetic oligo (dT) beads and then divided into short fragments by fragmentation buffer (Ambion, Austin, TX, United States). First-strand cDNA was synthesized with random hexamer primers, and then second-strand cDNA was synthesized using DNA polymerase I ( To facilitate assembly, the following criteria were used to filter low-quality reads: 1. Filter reads with adapter contamination. 2. Filter reads with unknown nucleotides >5%. 3. Filter reads in which more than 20% of bases showed a Q-value less than 20. High-quality clean reads was mapped to the reference genome to calculate the expression level in Fragments Per Kilobase of transcript per Million fragments mapped (FPKM).

Yeast Two-Hybrid Assays
Yeast two-hybrid (Y2H) assays were performed using the Matchmaker Gold Yeast Two Hybrid System (Clontech). The coding regions of the PaAQP genes (Supplementary Table 1) and predicted interacting proteins (Supplementary Table 2) were cloned into pGBKT7 and pGADT7 vectors to create different baits and prey, respectively. The different pairs of bait and prey constructs were then cotransformed into yeast strain Gold Y2H, and yeast cells were grown on plates lacking leucine and tryptophan (SD/-Leu/-Trp) for 3 days. Transformed colonies were plated onto quadruple dropout medium (SD/-Leu/-Trp/-His/-Ade) containing 100 µM aureobasidin A and 3 mg mL −1 X-α-Gal at 30 • C to test for interactions between PaAQPs and predicted interacting proteins.

Gene Expression Analysis and Quantitative Real-Time-PCR
To study the gene expression profiles of the PaAQPs in three periods of dormancy and germination in P. armeniaca 'Youyi, ' RNA sequencing (RNA-seq) data were collected from the transcriptome sequencing data from flower buds collected during the ED stage, SP, and GS of the of two Youyi plants grown under the same conditions.
The bud and stem samples at each stage were collected according to Song et al. (2017). Flower buds were checked and did not burst under forced conditions on the sampling day of PD or ED period, and could burst under forced conditions on the sampling day of SP or GS period. When flower buds are in the SP period, flower buds are not covered by scales, and petals are not visible. While flower buds are in the GS period, they are not covered by scales, and petals are visible. The RNA-seq data of P. armeniaca 'Youyi' in this study were deposited in the NCBI Sequence Read Archive (SRA) database under number SRS1042411. On the basis of fragments per kilobase of transcript per million mapped reads (FPKM) values, heatmaps and hierarchical clusters were analyzed with HemI 1.0 18 . Quantitative real-time (qRT)-PCR detection was chosen to analyze the expression of PaAQP genes in flower buds and stems of three P. armeniaca plants exhibiting similar growth during the PD stage, ED stage, SP, and GS. Total RNA was extracted from flower buds and stems using the RNeasy Plant Mini Kit (Qiagen, Valencia, CA, United States), and reverse transcription reactions were performed using 18 http://hemi.biocuckoo.org/down.php Superscript II reverse transcriptase (Invitrogen, Grand Island, NY, United States) according to the manufacturer's protocols. The primers used for qRT-PCR were designed using Primer Premier 5.0 software (Premier Biosoft Int., Palo Alto, CA, United States), and the specific primer sequences are listed in Supplementary Table 3. qRT-PCR was performed with an ABI Prism 7500 sequence detector (Applied Biosystems, Foster City, CA, United States) and the SYBR R Premix Ex TaqTM Kit (TaKaRa, Tokyo, Japan) to measure transcript levels. PaElf was used as a reference gene in qRT-PCR detection. qRT-PCR analyses were performed in triplicate for each biological sample, and quantitative results were analyzed according to the 2 − CT method (Livak and Schmittgen, 2001).

Coexpression Network Analysis
The coexpression data for PaAQPs were obtained from RNA-Seq data (SRA database number SRS1042411). For P. armeniaca genome-wide coexpression network construction, transcriptome data from the flower buds of two P. armeniaca plants were used. Genes with Pearson correlation coefficients >0.876 were selected, and 1,500 genes were screened preliminarily. Among the 1,500 genes coexpressed with PaAQPs, 14 genes related to cold resistance and 16 genes associated with cold stress were selected to build the coexpression network of PaAQPs. A network diagram was generated according to Cytoscape software (The Cytoscape 1.1 Core was downloaded from 19 ).

Yeast Transformants and Low-Temperature Treatment
The coding DNA sequences (CDSs) of PaPIP1-1, PaPIP1-3, PaPIP2-3, or PaTIP1-1 were inserted into the yeast expression vector pGAPZA. The recombinant vector pGAPZA:PaAQP was introduced into yeast strains. After identifying the positive colonies by PCR, the low-temperature treatment was conducted in solid and liquid Simmons Citrate-Ura (SC-U) media as described previously (Liu et al., 2015). Briefly, colonies of positive clones harboring the PaAQP gene or empty-vector (pGAPZA vector) were shaken and cultured at 200 rpm and 30 • C for 24 h to induce PaAQP protein expression. A 1-mL aliquot of yeast culture liquid (containing the pGAPZA empty or pGAPZA:PaAQP recombinant vectors) was added to the strains under the low-temperature treatment (at an OD600 value of 1.0), which was then cultured at −20 • C for 24 h. Afterward, the yeast liquid was diluted 1:10, 100, 1,000, and 10,000, followed by the sequential addition of 4 µL of the original yeast liquid and the diluted liquid to solid SC-U media (consisting of 2% galactose). After 48 h of culture at 30 • C, the growth rate of the transformed yeast cells was observed and recorded. In addition, 1 mL of the aliquot yeast (containing the pGAPZA empty vector or pGAPZA:PaAQP recombinant vector) liquid at an OD600 value of 1.0 was removed and added to 10 mL of liquid SC-U medium (containing 2% galactose). After 24 h of −20 • C treatment and 200 rpm at 30 • C, the OD600 value of the transformed yeast cells was measured. A control group that was not stressed was also included. 19 http://www.cytoscape.org/

Plant Transformation and Low-Temperature Resistance Analysis
The PaPIP1-3 and PaTIP1-1 CDSs were cloned into a pBI121 plant expression vector. The pBI121:PaPIP1-3 and pBI121:PaTIP1-1 recombinant plasmids were introduced into Agrobacterium tumefaciens (LBA4404). Transformation into Arabidopsis was carried out via the floral-dip method (Clough and Bent, 2010). To test the effect of cold stress on seedling growth, the seeds were allowed to germinate under normal growth conditions and then transferred to growth chambers at 22, 16, and 16 • C 16 h/4 • C 8 h. After 10 days of treatment, the physiological indexes of transgenic Arabidopsis and wild type plants were determined. Superoxide dismutase (SOD) activity (Wang et al., 2012a) and malondialdehyde (MDA) and proline levels were measured as previously described (Bates et al., 1973). Each sample comprised at least three individual plants, and all physiological analysis experiments were conducted in triplicate.

Statistical Analyses
The data concerning plant height, gene expression, OD600 values, etc., were assessed by one-way analysis of variance (ANOVA) followed by Tukey's HSD post hoc analysis. All the tests were conducted in triplicate. In the figures, the lowercase letters indicate statistical significance based on one-way ANOVA (P < 0.05). Statistical analyses were performed with SPSS 16.0 software (SPSS, Inc., Chicago, IL, United States).

Identification, Classification, and Properties of PaAQP Genes in P. armeniaca
Based on an extensive analysis of AQP-homologous sequences, a total of 33 putative full-length cDNA sequences of AQPs were identified from P. armeniaca. To investigate the phylogenetic relationship of AQP proteins among P. armeniaca, A. thaliana, and P. persica, a phylogenetic tree was constructed based on 33 protein sequences of AQP genes obtained from P. armeniaca, 35 sequences from A. thaliana, and 29 sequences from P. persica using the maximum likelihood method. The PaAQPs were classified into five subfamilies: NIPs (10 members), SIPs (4 members), TIPs (10 members), XIPs (2 members), and PIPs (7 members). The NIP subfamily included seven subgroups (one NIP1, one NIP2, two NIP3s, two NIP4s, two NIP5s, one NIP6, and one NIP7). The TIP subfamily was divided into five subgroups (three TIP1s, three TIP2s, two TIP3s, one TIP4, and one TIP5). The SIP subfamily was divided into two subgroups (three SIP1s and one SIP2s). The PIP subfamily was further clustered into three PIP1 and four PIP2 subgroups, and two members (XIP1s) belonged to the XIP subfamily ( Figure 1A and Supplementary Table 1). All of these genes were grouped with A. thaliana and P. persica AQPs, confirming that they belonged to the AQP family.
Synteny analysis was carried out using MCScanX software, and syntenic relationships and whole genome sequences to visualize the locations of orthologous genes were presented using TBtools software. Microsynteny analysis was performed across the P. armeniaca, A. thaliana, and P. persica genomes. Between P. armeniaca and A. thaliana, 14 collinear blocks were identified, while 22 orthologous gene pairs were found ( Figure 1B). Between P. armeniaca and P. persica, 30 collinear blocks were identified, with 41 orthologous gene pairs ( Figure 1B). These results suggest that P. armeniaca and P. persica are more closely related than P. armeniaca and A. thaliana.
Bioinformatics analysis, CDS and amino acid sequence analysis, and ExPaSy database analysis showed that the Mw of the PaAQPs ranged from 21.94 to 35.37 kDa and that their pIs ranged from 5.32 to 10.00 ( Table 1). The SIPs and TIPs were smaller (<28 kDa) than the NIPs, XIPs, and PIPs. The TIPs were acidic (except for PaTIP2-3), while the majority of the other subfamilies were alkaline. Most of the AQPs were predicted to have TMHs, while PaTIP2-2 and PaTIP2-3 had only five, PaNIP3-1 had seven, and PaNIP5-1b had eight TMHs ( Table 1). The GRAVY scores of the PaAQP proteins ranged from 0.291 to 0.990, and among Frontiers in Plant Science | www.frontiersin.org FIGURE 1 | Phylogenetic relationships and microsynteny analyses of AQP genes among the P. armeniaca, A. thaliana, and P. persica genomes. (A) Deduced amino acid sequences were aligned using default parameters in ClustalW, and the phylogenetic tree was constructed via the maximum likely hood method with 1,000 bootstrap replicates with MEGA 7.0 software. Red, green, and blue circles represent the A. thaliana, P. armeniaca, and P. persica AQP gene family members. Yellow, blue, red, orange, and green arches represent the SIP, XIP, NIP, PIP, and TIP AQP gene family subgroups of A. thaliana, P. armeniaca, and P. persica. The scale bar indicates the distance calculated by way of multiple alignment. (B) The chromosome numbers of all three species are specified by different colors: blue, green, and yellow represent the P. armeniaca, A. thaliana, and P. persica chromosomes, respectively. The chromosome number is indicated on the inside with the chromosome sequence lengths in megabases. Gene pairs with syntenic relationships are linked by blue and green lines representing the microsyntenic regions between the P. armeniaca and A. thaliana chromosomes and the P. armeniaca and P. persica chromosomes, respectively. the subfamilies, the PIPs presented the lowest average GRAVY value (0.409) ( Table 1).
Gene duplication has long been considered one of the main forces in the evolution and expansion of a gene family (Crow and Wagner, 2006). Therefore, we investigated the different types of gene duplications in the PaAQP gene family. We found that 6 AQP genes clustered in 3 tandem repeat event regions in the P. armeniaca genome (Supplementary Figure 1). There were three pairs of tandem duplicated genes (PaNIP5-1a and PaNIP5-1b, PaTIP2-2 and PaTIP2-3, PaXIP1-1, and PaXIP1-2) within the 1 Mb chromosomal region of chromosomes 3, 5, and 8. In addition, 6 pairs of WGD/segmental duplication genes (PaTIP4-1 and PaTIP2-1, PaTIP1-2 and PaTIP1-3, PaSIP1-1, and PaSIP1-2, triangle WGD/segmental duplication among PaPIP2-1, PaPIP2-2, and PaPIP2-3) were found within the PaAQP family (Supplementary Figure 1). Notably, all the identified segmental duplication gene pairs were distributed on different chromosomes in the P. armeniaca genome. However, the distribution of segmental replication genes on the chromosomes was uneven; there were three genes on chromosome 6, two on chromosomes 2 and 7, and only one on chromosome 1. Moreover, we also found that PaPIP2-1 was the common duplicated gene between PaPIP2-2 and PaPIP2-3. This particular gene replication event involved three genes: PaPIP2-1, which was located on chromosome 2, and PaPIP2-2 and PaPIP2-3 on chromosome 6. Additionally, the analysis revealed one particular WGD/segmental duplication event between chromosomes 2 and 6 involving four genes. PaPIP2-1 combined with PaTIP1-3 was reverse duplicated to PaPIP2-2 combined with PaTIP1-2.
Exon-intron organization plays key roles in phylogenetic relationships and plant evolution Jue et al., 2015). As shown in Supplementary Figure 2, the distribution of introns in the PaAQP genes ranged from one to five introns per gene. All PaPIPs presented four exons: PaXIP1-1 had two exons, and PaXIP1-2 had three exons. The majority of PaSIPs exhibited three exons; an exception involved PaSIP1-3, which had only one exon and was therefore intronless. Most PaTIPs had three exons, while PaTIP1-1 had two and PaTIP2-3 had five. A relatively diverse exon-intron pattern was detected for the PaNIP subfamily. The majority of PaNIPs had five exons, with the exceptions of PaNIP3-2, which had six exons, PaNIP5-1a, which had four exons, and PaNIP3-1, which had three exons. Exonintron organization analysis revealed a conserved pattern of gene structure within the subfamilies of PaAQPs. The conserved exon-intron structure supports the phylogenetic relationships of P. armeniaca (Supplementary Figure 2).

Characterization of Asn-Pro-Ala Motifs, Transmembrane Domains, and Conserved Motifs of PaAQPs
To better understand the possible physiological mechanisms of the PaAQPs, all 33 sequences were aligned, and the conserved residues (NPA motifs, ar/R selectivity filter, and Froger's position) were analyzed (Figure 2, Table 2, and Supplementary  Figures 3-6). The P. armeniaca PaAQPs displayed differences in NPA motifs and residues at the ar/R selectivity filters and Froger's positions compared with those of other plant species. Most of the AQPs showed dual NPA motifs, with the exception of PaSIP2-1, which was found to harbor a single NPA motif. The majority of members from the PIP and TIP subfamilies contained a typical NPA motif, as observed for the Beta vulgaris counterpart (Kong et al., 2017), with the exception of PaTIP2-3, which showed an alanine (A)-to-serine (S) substitution in the first NPA motif and asparagine (N)-to-A and A-to-lysine (K) substitutions in the secondary NPA motif. The LPK (PaTIP2-3) motifs detected in P. armeniaca have not been reported in any other plant species. These changes may alter the substrate specificity of the AQPs in P. armeniaca. In the NIP subfamily, some members of the NIPs, such as PaNIP5-1a, PaNIP3-2, PaNIP5-1b, and PaNIP6-1, showed replacement of alanine (A) at the third residue of the first NPA motif by serine (S) and replacement of other alanine FIGURE 2 | Multiple sequence alignment of the deduced amino acid sequences of PaPIPs. Amino acids at NPA domains, ar/R selectivity filters, and Froger's residues were identified in seven PIP family members (3 PIP1s and 4 PIP2s) in P. armeniaca.
A) residues by isoleucine (I) or valine (V). In the XIP subfamily, the first NPA motif showed an N-to-S substitution in PaXIP1-1 and a valine (V)-to-A substitution in PaXIP1-2, while the second NPA motif was conserved. All PIP subfamily members showed conserved ar/R filter residues, with phenylalanine (F) at H2, histidine (H) at H5, threonine (T) at LE1, and arginine (R) at LE2, which is typical of water-transporting AQPs. They also showed Q/E-S-A-F-W residues at Froger's positions, as observed in species such as flax (Linum usitatissimum) (Shivaraj et al., 2017) and chickpea (Cicer arietinum L.) (Deokar and Bunyamin, 2016). In the TIP subfamily, the H2 position of the ar/R filter contained histidine (H), and the H5 position contained isoleucine (I), except in PaTIP1-1, where this position was occupied by valine (V), and PaTIP2-3, which did not exhibit a residue at this position. The LE1 and LE2 positions were found to be specific to each subgroup of PaTIP. The PaTIP1 subgroup was characterized by alanine (LE1) and valine (LE2) levels. The PaTIP2 subgroup was characterized by glycine (LE1) and arginine (LE2), except for PaTIP2-2, which contained glycine (LE2), and PaTIP2-3, which contained alanine (LE1) and leucine (LE2). The PaTIP3 and PaTIP4 subgroups were characterized by alanine (LE1) and arginine (LE2), with the exception of PaTIP3-1, which contained glycine (LE1) and cystine (LE2). PaTIP5 was also characterized by glycine (LE1) and cystine (LE2) and showed a conserved relationship with chickpea (Deokar and Bunyamin, 2016) and flax (Shivaraj et al., 2017). Among the NIPs, the ar/R selectivity filter and Froger's positions presented multiple types. PaNIP1-1, PaNIP3-1, PaNIP4-1, and PaNIP4-2 showed W-R To further analyze the conservation, diversity, and relationship in each subfamily, motif scan analysis was conducted via the MEME program. Seventeen conserved motifs were identified in the PaAQP family, and most of the members in the same subfamily shared a similar number of motifs (Supplementary Figure 7). For example, all the members of PaPIPs had thirteen motifs, except for PaPIP1-2. All the members of PaTIP1s had eleven motifs, all the members of PaSIP1s had nine motifs, and all the members of PaNIP4s had eleven motifs. The results also showed that PaAQPs had similar functional motifs (motif 2) consisting of highly conserved regions (channel proteins of the major intrinsic protein (MIP) superfamily) with motif 4 from B. vulgaris BvAQP proteins. These results implied that PaAQPs may have a similar function to AQPs from other plants.

Molecular Modeling of PaAQP
3D proteins of plant AQPs can be effectively used for understanding the substrate specificity and solute permeation rate of AQP genes, which will help to improve the understanding of the cold resistance of PaAQPs (Groot and Grubmüller, 2005). The structural properties of all 33 PaAQPs were displayed in homology-based tertiary (3D) protein models, which were predicted via the Phyre2 server display, and the results are shown in Supplementary Table 5. All 3D protein models were constructed with 100% confidence, and the residue coverage varied from 68 to 98%. The PaAQP 3D protein model contained a conserved hourglass-like structure consisting of an a-helical bundle forming five to eight TM helices (H1 to H8) and two short helices (HE and HB). Loops HE and HB were located close to each in the center of the membrane. These predicted 3D models showed the boron (B) and silicon (Si) substrate selectivity of the PaAQPs and provided an important basis for the protein functional analysis of the AQP genes of P. armeniaca (Deshmukh et al., 2013;Tombuloglu et al., 2016).

Protein-Protein Interaction Network of PaAQPs
To further reveal the function of PaAQPs during the interaction with other proteins, PaTIP1-1, PaPIP1-3, PaSIP1-3, and PaNIP6-1, which showed high expression in each subfamily (Figure 3), were used to construct a protein-protein interaction network. As shown in Figures 4A-D and Supplementary Table 6, none of these tested PaAQPs shared the same interacting proteins. In addition, various proteins specifically interacted with TIP separately, such as zinc finger proteins, calciumdependent protein kinases, and NAC domain-containing proteins. In particular, the interacting zinc finger proteins and NAC transcription factors can enhance cold tolerance in transgenic plants . Interestingly, TIP2 FIGURE 3 | Expression profiles of PaAQPs at different developmental stages based on RNA-seq data. ED, SP, and GS indicate that the tested materials were collected during ecological dormancy, the sprouting period, and germination stage of P. armeniaca, respectively. The number in the middle of each box represents the FPKM. could interact with TIP and form homodimers. The proteins that interact with PIP include AGD2-like defense response protein, methylenetetrahydrofolate reductase family proteins, and Delta 1-pyrroline 1-5-carboxylate reductase. Specifically, SIP and NIP members could interact with each other. For example, NIP6-1 specifically interacted with SIP1-3, SIP2-1, and SIP1A, while four NIP members (NIP1-2, NIP2-1, NIP4-1, and NIP6-1) and three PIP members (PIP1B, PIP1B, and PIP2-5) could interact with SIP1A. Notably, NIP6-1 specifically interacted with WRKY6, which plays important roles in cold resistance, stress tolerance, growth and development . Y2H assays showed that yeast cells cotransformed with constructs PaNIP6-1-pGBKT7 and WRKY6-pGADT7, PaTIP1-1-pGBKT7 and TSPO-pGADT7, and PaPIP2-1-pGBKT7 and ATHATPLC1G-pGADT7 could support growth in the selection medium and b-galactosidase staining ( Figure 4E). Thus, Y2H assays revealed the protein interaction between PaNIP6-1 and WRKY6, PaTIP1-1 and TSPO, and PaPIP2-1 and ATHATPLC1G. In particular, the A. thaliana gene encoding a phosphatidylinositol-specific phospholipase C AtPLC1 (ATHATPLC1G) was induced by low-temperature, dehydration, and salt stress, which might contribute to the adaptation of the plant to cold stresses (Hirayama et al., 1995).

Expression Profiling of PaAQPs in the Dormancy and Sprouting Stages of P. armeniaca
To gain more information about the potential role of PaAQP proteins in P. armeniaca, we analyzed the expression of PaAQPs in the ED, SP, and GS by transcriptome sequencing. The RNAseq data can be found in the NCBI SRA database under number SRS1042411. The expression levels of 33 PaAQP genes were clustered in a heatmap (Figure 3). We observed the expression of 30 PaAQPs (90.90%) in at least one of the developmental stages, and all of these PaAQPs showed differential expression in three stages of dormancy and sprouting periods. Among these genes, 23 PaAQP genes (75%) were expressed in all three periods of dormancy and germination, among which PaPIP2-1, PaPIP1-3, and PaTIP1-1 exhibited the greatest expression levels in different periods. In addition, PaTIP1-2 showed greater expression in the ED period than in the SP and GS periods, which may suggest that the P. armeniaca plant was in the ED period. Similarly, PaNIP7-1 showed greater expression in the SP period than in the ED and GS periods and may be useful as a marker gene of entry into the SP. PaTIP1-3 and PaNIP2-1 exhibited greater expression levels in the GS than in the ED and SP, so these two genes might be useful for suggesting that the P. armeniaca plant is entering the GS. Interestingly, PaNIP5-1b, PaTIP2-3, and PaTIP2-2 were not expressed in any of the three periods of dormancy and germination. In addition to these three genes, PaNIP3-1, PaNIP4-1, PaNIP4-2, PaNIP5-1a, PaTIP3-1, and PaTIP5-1 were not expressed in the ED period of flower buds; and PaNIP5-1a, PaNIP4-1, and PaTIP3-1 were not expressed in the SP of flower buds.

Coexpression Network of PaAQPs in P. armeniaca
To elucidate the molecular mechanism of PaAQPs in the dormancy and sprouting stages of P. armeniaca, a PaAQP coexpression network was constructed (Figure 6A and FIGURE 5 | Changes in the transcript levels of 20 selected genes in different tissues during the dormancy and sprouting stages of P. armeniaca. The stem and flower buds were collected during the PD stage, ED stage, SP stage, and GS stage of P. armeniaca plants. All expression levels of PaAQP genes were normalized to the expression levels of PaElf (elongation factor-1α). The 2 − CT method was used to analyze the relative gene expression. The data are means of three replicates ± SDs. The lowercase letters indicate statistical significance based on one-way ANOVA with Tukey's HSD post hoc analysis.
FIGURE 6 | PaAQP coexpression network. (A) Among the 1,500 genes coexpressed with PaAQPs, 14 (green nodes) are involved in cold resistance based on previous research. In addition, the anticold ability of 16 genes (pink nodes) is not very clear. The red node is the PaAQP genes. Cyan lines represent PaAQPs coexpressed with 14 cold resistance genes, and red lines represent PaAQPs coexpressed with 16 cold stress-associated genes. (B) Transcript levels of 14 genes involved in cold resistance in stem and flower buds during the PD stage, ED stage, SP stage, and GS stage of P. armeniaca. All expression levels of the PaAQP genes were normalized to those of PaElf. Each bar represents the mean ± SD of three technical replicates. The different letters represent significant differences at P < 0.05 (one-way ANOVA).

Supplementary Table 7).
A total of 1,500 genes were identified as coexpressed with PaAQPs. Among these, 14 genes (0.93%) were associated with cold resistance. For example, C3HC4-type RING zinc finger proteins are involved in cold stress . MYB transgenic rice enhances tolerance to chilling stress and mediates ectopic expression of stress genes (Ma et al., 2009). bZIP transcription factors regulate freezing tolerance in Arabidopsis (Yao et al., 2020). Additionally, overexpression of the NAC transcription factor enhances the cold tolerance of Medicago truncatula (Qu et al., 2016). Serine/threonine protein kinases are involved in the cold stress response in the cyanobacterium Synechocystis (Zorina et al., 2014). In our research, zinc finger proteins, MYBs, bZIPs, NACs, and serine/threonine protein kinases, among others, were highly coexpressed with PaAQPs in different tissues during the dormancy and sprouting stages of P. armeniaca and may be involved in cold resistance in P. armeniaca.
To validate the network constructed by the core PaAQPs, changes in the expression of the 14 genes related to cold resistance in the PaAQP coexpression network were detected in the PD, ED, SP, and GS of flower buds and stems of P. armeniaca plants ( Figure 6B). Among the 14 pairs of coexpressed genes, the expression trends of PaSIP1-1 and PaZFC3HC4, PaSIP1-1 and PaEREBP were consistent in the four periods (PD, ED, SP, and GS) of flower buds of P. armeniaca, whereas those of PaTIP4-1 and PaSGBP were rarely consistent. Similarly, the expression was similar between PaNIP6-1 and PaUbI and between PaNIP6-1 and PaWRKY. Additionally, three pairs exhibited no differences in expression patterns: PaXIP1-2 and PaMYB, PaTIP1-3 and PaPRE, PaPIP2-4 and PaLEA.
Subcellular Localization of PaPIP1-1, PaPIP2-3, PaSIP1-3, PaXIP1-2, PaNIP6-1, and PaTIP1-1 Proteins WoLF PSORT and Plant-mPLoc prediction were used to predict the subcellular localization of the PaAQPs. Almost all members of the PaPIP subfamilies were localized to the plasma membrane according to WoLF PSORT and Plant-mPLoc prediction (Table 1). However, the subcellular localizations predicted by Plant-mPLoc indicated that PaPIP2-2 was localized to the cell wall. The PaTIPs were predicted to localize to the vacuoles, and PaNIPs (except for PaNIP5-1a),PaSIPs, and PaXIPs were predicted to localize to the plasma membrane according to Plant-mPLoc prediction. However, the results predicted by WoLF PSORT were diverse and included localization to the nucleus, Golgi, cytosol, mitochondrion, chloroplast, and endoplasmic reticulum ( Table 1). Genes with high expression in different periods of low temperature may play important roles in cold resistance during the dormancy and sprouting stages of P. armeniaca. To identify and gain insight into the subcellular localization of these highly expressed PaAQP genes in plant cells, chimeric EGFP-PaPIP1-1, EGFP-PaPIP2-3, EGFP-PaSIP1-3, EGFP-PaXIP1-2, EGFP-PaNIP6-1, and EGFP-PaTIP1-1 cDNAs (representing genes with the greatest expression in each subgroup in flower buds in the dormancy and germination of P. armeniaca) were constructed and inserted downstream of the CaMV35S promoter. Transient expression analyses were performed using A. thaliana leaf mesophyll cells.
Analyses of fusion protein localization by confocal microscopy showed that the green fluorescence of these proteins was confined to the plasma membrane (Figure 7). This membrane localization of the PaPIP1-1, PaPIP2-3, and PaSIP1-3 proteins was consistent with the location of PIP and SIP gene expression in other plants (Zou et al., 2015;Deokar and Bunyamin, 2016).
The PaPIP1-3 and PaTIP1-1 Genes Conferred Low-Temperature Resistance to Yeast To analyze the function of PaAQP genes involved in lowtemperature tolerance, four genes that showed high expression in flower buds and stems during the dormancy and sprouting stages of P. armeniaca were expressed in yeast strain GS115 to check their low-temperature stress resistance. The yeast strains transformed with the pGAPZA empty vector or various pGAPZA:PaAQP recombinant expression vectors were incubated in solid and liquid SC-U medium (2% galactose) for low-temperature stress treatment. The growth of yeast cells transformed with pGAPZA:PaPIP1-3 and pGAPZA:PaTIP1-1 on solid SC-U media consisting of 2 M sorbitol after −20 • C low-temperature treatment grew much more than did the cells transformed with empty vector. Moreover, the colonies of yeast cells harboring the PaPIP1-3 and PaTIP1-1 genes were larger than those harboring empty vectors strains ( Figure 8A). In 2 M liquid SC-U sorbitol media, the mean OD value of the empty-vector strains was 2.33, whereas the positive strains harboring PaPIP1-3 and PaTIP1-1 had OD values of 2.73 and 2.55, respectively. After the −20 • C low-temperature treatment, the mean OD values were 1.99 vs. 2.54 and 2.48 ( Figure 8B). These results suggested that the PaPIP1-3 and PaTIP1-1 genes improved the low-temperature stress resistance of yeast cells.

PaTIP1-1 in A. thaliana Conferred Cold Tolerance
To study the function of the PaAQP genes in cold stress, two genes (PaPIP1-3 and PaTIP1-1) exhibiting cold stress resistance in yeast were introduced individually into A. thaliana plants. Normally growing wild-type Arabidopsis and PaPIP1-3 or PaTIP1-1 gene-carrying Arabidopsis were transferred to growth chambers at 22, 16, or 16 • C 16 h/4 • C 8 h. After 10 days of treatment, the wild-type plants showed significantly high sensitivity to cold stress with yellow leaves. The growth of PaPIP1-3-or PaTIP1-1-overexpressing (OE) Arabidopsis lines (PIPOE-1, PIPOE-2, TIPOE-1, and TIPOE-2) was better than that of the wild-type plants, demonstrating good cold resistance ( Figure 9A). After low temperature treatment, two independent PaPIP1-3 transgenic lines (PIPOE-1 and PIPOE-2) and two independent PaTIP1-1 transgenic lines (TIPOE-1 and TIPOE-2), which with higher relative expression level than untreated transgenic lines, were identified by qRT-PCR ( Figure 9B). We further examined the activities of SOD, the major ROS antioxidant enzyme in plants. The PIPOE-1 and PIPOE-2 lines and the TIPOE-1 and TIPOE-2 lines had significantly greater SOD activity than the wild-type line ( Figure 9C). The superior performance of the PaPIP1-3-and PaTIP1-1-transformed plants was further validated by analyzing the proline and MDA contents subsequent to cold stress. The proline level was significantly greater in the transgenic leaves than in the wild-type leaves under the stress conditions, indicating superior biochemical capabilities of the transgenic plants ( Figure 9D). Similarly, the MDA content was greatest in wild-type leaves of cold-stressed plants, showing that the greatest membrane damage had occurred in untransformed plants (Figure 9E). This research indicated that the target genes PaPIP1-3 or PaTIP1-1 were directly related to cold resistance, and the increase in PaPIP1-3 or PaTIP1-1 expression in PaPIP1-3-or PaTIP1-1-transformed plants could improve the physiological condition and cold resistance of transgenic plants.

DISCUSSION
Aquaporins are varied and are widely distributed and functionally diverse in plant genomes. These proteins play important roles in the transport of water and substrates in plants, in regulating plant development and growth, and in responses to drought, cold, or salt stresses (Prado and Maurel, 2013;Kong et al., 2017). A majority of the genome-wide analyses and functional studies of AQP genes conducted to date have focused on model plants such as Z. mays, A. thaliana, O. sativa, N. tabacum, and Populus (Chaumont et al., 2001;Gupta and Sankararamakrishnan, 2009;Nguyen et al., 2013). Very little is known about the AQPs of economic forest species such as P. armeniaca. To determine the evolutionary relationships and functional mechanisms of PaAQPs, a phylogenetic tree was built using AQPs from A. thaliana (in which a full set of AtAQP genes is well known) and Prunus persica (Deshmukh et al., 2015), which are closely related species of P. armeniaca. The number of AQPs identified in this study was the same as that reported in Brachypodium distachyon (Deshmukh et al., 2015), being slightly greater than the 30 AQP genes reported in Elaeis guineensis (Deshmukh et al., 2015) and Vitis vinifera (Deshmukh et al., 2015) and lower than the 33 genes found in rice (Sakurai et al., 2005), 35 genes found in A. thaliana (Zhang et al., 2013), and 37 genes found in Ricinus communis (Deshmukh et al., 2015). The number of P. armeniaca AQPs in the five subfamilies was also similar to their distribution in Cajanus cajan (Deshmukh et al., 2015). However, the numbers of PIPs and TIPs varied in the two genomes compared with those of the other subfamilies. Twelve and thirteen members of PIPs and TIPs have been reported in C. cajan (Deshmukh et al., 2015), while seven SIPs and ten TIPs have been detected in P. armeniaca. PaTIPs are divided into five subgroups, TIP1, TIP2, TIP3, TIP4, and TIP5, and this classification is consistent with those previously reported in species such as barley, sugar beet, rice, and maize (Tombuloglu et al., 2016;Kong et al., 2017).
of the PaAQP genes were closely related to the AQP genes in P. persica, consistent with the finding that P. armeniaca and P. persica diverged from a recent common ancestor. It is well known that gene duplication mechanisms (WGD/segmental duplication, tandem duplication) have a significant role in biological evolution (Panchy et al., 2016). In the present study, we observed that orthologous genes developed through segmental duplication and tandem duplication. Both styles of duplication have played a significant role in the expansion of P. armeniaca AQP family genes. The identification of orthologous PaAQP genes will further benefit the evolutionary history of this gene family.
Asn-Pro-Ala motifs, the ar/R selectivity filter, and Froger's positions are considered to participate in transport activity and substrate selection. PaAQP functions can be speculated based on the amino acid residues of the proteins compared with the AQPs of other plants. The PaPIPs showed typical NPA motifs, F-H-T-R residues in the ar/R selectivity filter, and Q/E-S-A-F-W at Froger's positions, which are highly conserved in the PIPs of other plants, such as flax (L. usitatissimum) (Shivaraj et al., 2017), chickpea (C. arietinum L.) (Deokar and Bunyamin, 2016) and H. brasiliensis (Murata et al., 2000). This composition of PIPs is considered to play a vital role in facilitating CO 2 diffusion, regulating leaf and root hydraulics, and influencing plant photosynthesis (Gupta and Sankararamakrishnan, 2009). Therefore, the homologous PIPs of P. armeniaca may play similar roles in regulating plant growth and development, such as in photosynthetic physiology, plant hydraulics, and water and solute transport. Plant TIPs have been suggested to function as water transporters as efficiently as PIPs (Murata et al., 2000). For example, PaTIP1-1, BvTIP1-2, and BvTIP1-3 contain dual conserved NPA motifs, H-I/V-A-V at the ar/R selectivity filter and T-S-A-Y-W at Froger's positions, which is the same composition observed for the Citrus sinensis AQPs CsTIP1s and CsTIP3s. Additionally, these proteins are predicted to function in the transport of urea and H 2 O 2 (Hove and Bhave, 2011). Compared with other TIPs, PaTIP1-2, PaTIP1-3, and PaTIP2-1 showed H-I-A-V or H-I-G-R at the ar/R selectivity filter and at Froger's positions and harbored T-S-A-Y-W residues at the P1-P5 positions; this structure type is considered to transport substances such as NH 4 , urea, and H 2 O 2 (Hove and Bhave, 2011). These analyses suggest that PaTIPs play an important role in the transport of a wide range of molecules.
WoLF PSORT predicted that all the PaPIP members were localized in the plasma membrane, and plant-mPLoc prediction indicated that the PaNIPs (except for PaNIP5-1a) were also localized in the plasma membranes. These results were consistent with experimental findings for PIPs and NIPs of banana, flax, and other plant species (Maurel et al., 2008;Shivaraj et al., 2017). However, some tobacco NtPIP members show dual localization in the plasma membrane and inner chloroplast membrane (Uehlein and Kaldenhoff, 2008), and different localizations of a maize PIP gene (ZmPIP1-2) have been found in the endoplasmic reticulum and plasma membrane (Chaumont et al., 2000). NIPs also exhibit different subcellular localizations in the peribacteroid membrane, endoplasmic reticulum, or plasma membrane in other plants (Takano et al., 2006;Deokar and Bunyamin, 2016). According to our predictions, both PaNIP1-1 and PaNIP4-2 were localized to the plasma membrane and vacuoles. The reasons for the different localizations of these genes are not clear. The majority of SIPs have been shown to localize to the vacuoles and plasma membrane; however, grapevine and A. thaliana SIPs are predicted to be localized to the endoplasmic reticulum Noronha et al., 2014;Kong et al., 2017). In addition to differences between species and protein sequences, many experiments have shown that external environmental conditions such as drought and salinity, but not cold, influence the characteristics of the subcellular localization of PaAQPs (Boursiac et al., 2008;Deokar and Bunyamin, 2016). These results indicate that the subcellular localization of plant AQPs is complex and regulated by environmental conditions. Aquaporin genes are expressed in different plant tissues and organs, and their spatiotemporally specific expression is highly correlated with different developmental stages (Alexandersson et al., 2005;Giorgio et al., 2016). The expression patterns of AQPs in various tissues as well as under different stresses may provide an efficient basis for identifying their molecular functions. qRT-PCR results showed that PaAQPs were expressed in all tissues tested, but their transcript levels largely differed in various organs (Figure 5), suggesting that the functions of PaAQPs might vary. For example, the expression levels during the endodormancy stage (PaPIP2-3, PaPIP2-2, and PaSIP1-3 are high) indicated that the effects of growth inhibition were controlled by other internal factors resulting from insufficient low-temperature accumulation. High expression levels during the ecodormancy stage (PaTIP1-2 and PaTIP2-1) could reflect the activation of growth through release from endogenous inhibitors but maintained suppression due to unsuitable conditions (insufficient high temperature accumulation). Expression during SP (PaSIP1-2, PaNIP3-2, and PaNIP6-1 are high) and GS (PaPIP1-1, PaPIP2-4, PaSIP1-1, PaSIP2-1, PaXIP1-2, PaNIP2-1, PaNIP1-1, PaNIP3-1, PaTIP1-1, PaTIP1-3, and PaTIP4-1 are high) showed morphological growth resumption, not (endo)dormancy release. Moreover, PaPIP1-1, PaPIP2-2, PaPIP2-3, and PaSIP1-1 were highly expressed in stems, where they may contribute to the storage of energy and materials needed for growth, development, and germination.
In the PaAQP coexpression network, similar expression patterns of gene pairs (such as PaSIP1-1 and PaZFC3HC4, PaSIP1-1 and PaEREBP, PaNIP6-1 and PaUbI and PaNIP6-1 and PaWRKY) during the dormancy and sprouting stages suggested that PaAQPs cooperated with various vital cold stressor flower bud dormancy/germination-related genes to regulate cold tolerance, bud dormancy, and bud break. For example, a novel zinc-finger protein (homologous gene of PaZFC3HC4) mediates ABA-regulated seed dormancy in Arabidopsis (He and Gan, 2004), and an ubiquitin-conjugating enzyme (homologous gene of PaUbI) has been found to be involved in tuber dormancy breaking (Liu et al., 2012). A total of 68 AP2/ERF genes were first identified in dormant Chinese cherry flower buds. Some PpcAP2/ERF TFs (homologous genes of PaEREBP) participate in and influence the dormancy transition through their responses to low temperature (Zhu et al., 2021). Chen et al. (2016) found that the expression levels of six WRKY genes increased during endodormancy and decreased during ecodormancy, indicating that these six WRKY genes may play a role in dormancy in peach. However, several coexpression pairs, such as PaSIP1-3 and PabZIP and PaTIP2-1 and PaNAC, exhibited different expression regularities, implying that these gene pairs may have different roles and functions in different periods during the dormancy and sprouting stages of P. armeniaca.
Aquaporins regulate long-distance water transport between cells and can mediate the osmotic pressure of the cytoplasm to appropriately regulate plant growth and development, responses to environmental stress, and stomatal movement, among others. PaTIPs can mediate rapid water transport between cells, remove accumulated H 2 O 2 , and improve water transport activity, which is necessary to end dormancy (Erez et al., 1998;Yooyongwech et al., 2008). In our study, several PaTIPs, such as PaTIP1-1, PaTIP1-2, and PaTIP4-1, showed greater transcript levels in all four dormancy and sprouting stages of P. armeniaca. PaTIP1-2 showed greater expression in the ED period than SP and GS, which may represent a response to cold stress to aid in the removal of the accumulated induced H 2 O 2 , thereby promoting the end of dormancy. In particular, overexpression of PaTIP1-1 increased cold tolerance and accumulation in yeast (Figure 8), which may help to improve water transport activity and promote normal growth under low-temperature stress. Previous research has indicated that AtPIP2-5 and AtPIP2-6 are upregulated by cold stress in A. thaliana plants (Alexandersson et al., 2005). Transgenic banana plants overexpressing MusaPIP1-2 show a high tolerance to cold conditions (Shekhawat and Ganapathi, 2013). A. thaliana plants overexpressing AtPIP1-4 or AtPIP2-5 display greater cold tolerance than wildtype plants (Afzal et al., 2016). In our study, PaPIP1-2 and PaPIP2-4 exhibited high expression in the dormancy and sprouting stages in the flower buds of P. armeniaca. Interestingly, yeast cells expressing PaPIP1-3 and PaTIP1-1 exhibited markedly enhanced growth compared with that of yeast cells transformed with the empty vector on SG-U medium with low-temperature treatment (Figure 8), and PaPIP1-3 or PaTIP1-1 overexpressed in Arabidopsis showed enhanced cold tolerance by improving antioxidative enzyme activities (Figure 9). These results may suggest that PaPIPs present stronger cold resistance and play a vital role in the maintenance of normal growth and development under cold conditions.

CONCLUSION
Overall, this work presents the first genome-wide study of the P. armeniaca AQP gene family, in which 33 non-redundant PaAQP genes were identified in the P. armeniaca genome and phylogenetically clustered into five distinct subfamilies. Synteny analysis identified 14 collinear blocks containing AQP genes between P. armeniaca and A. thaliana and 30 collinear blocks between P. armeniaca and P. persica. The chromosome localization, gene structure, protein characteristics, conserved functional motifs, and homologybased 3D models of the PaAQP genes were further examined. Cis-motif analysis of the 2.0-kb sequences upstream of the AQP genes indicated the existence of several light-, hormone-, and stress-responsive and developmental stagespecific elements in the PaAQPs. The subcellular localization of PaPIP1-1, PaPIP2-3, PaSIP1-3, PaXIP1-2, PaNIP6-1, and PaTIP1-1 was examined, and the green fluorescence of the recombinant proteins was observed in the plasma membrane of isolated protoplasts. In silico analysis of proteinprotein interactions suggested the association of PaAQP proteins with specific cold tolerance and development pathways. Gene expression profiling of PaAQP showed that PaPIP1-3, PaPIP2-1, and PaTIP1-1 were highly expressed in flower buds during the dormancy and sprouting stages of P. armeniaca. Significantly, PaPIP1-3 and PaTIP1-1 conferred low-temperature stress resistance in a yeast expression system. Overexpression of PaPIP1-3 and PaTIP1-1 in Arabidopsis plants increased antioxidative enzyme activities and improved transgenic plant tolerance to cold stress. This study provides an important resource for the future functional analysis of PaAQPs and a theoretical basis for the breeding of new varieties via the genetic engineering of P. armeniaca.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
SL and T-NW conceived and designed the research. SL, YX, and SW performed all the experiments. YZ, GZ, XZ, and JL contributed reagents, materials, and analysis tools. JZ, XG, LW, and WM analyzed the data. SL wrote the manuscript. All authors read and approved the final manuscript. , and other cis-acting elements in the promoter region of P. armeniaca PaAQP genes from each AQP family. The cis-acting elements were determined in the 2-kb promoter region upstream of the translation initiation site using the PLACE database.
Supplementary Figure 9 | qRT-PCR analysis of the relative expression of 20 genes during the dormancy and sprouting stages of P. armeniaca.
Supplementary Table 1 | Complete coding sequences and corresponding amino acid sequences of the AQP genes of P. armeniaca.