Identiﬁcation of Aquaporin Gene Family in Response to Natural Cold Stress in Ligustrum × vicaryi Rehd.

: Plants are susceptible to a variety of abiotic stresses during the growing period, among which low temperature is one of the more frequent stress factors. Maintaining water balance under cold stress is a difﬁcult and critical challenge for plants. Studies have shown that aquaporins located on the cytomembrane play an important role in controlling water homeostasis under cold stress, and are involved in the tolerance mechanism of plant cells to cold stress. In addition, the aquaporin gene family is closely related to the cold resistance of plants. As a major greening tree species in urban landscaping, Ligustrum × vicaryi Rehd. is more likely to be harmed by low temperature after a harsh winter and a spring with ﬂuctuating temperatures. Screening the target aquaporin genes of Ligustrum × vicaryi responding to cold resistance under natural cold stress will provide a scientiﬁc theoretical basis for cold resistance breeding of Ligustrum × vicaryi . In this study, the genome-wide identiﬁcation of the aquaporin gene family was performed at four different overwintering periods in September, November, January and April, and ﬁnally, 58 candidate Ligustrum × vicaryi aquaporin (LvAQP) genes were identiﬁed. The phylogenetic analysis revealed four subfamilies of the LvAQP gene family: 32 PIPs, 11 TIPs, 11 NIPs and 4 SIPs. The number of genes in PIPs subfamily was more than that in other plants. Through the analysis of aquaporin genes related to cold stress in other plants and LvAQP gene expression patterns identiﬁed 20 LvAQP genes in response to cold stress, and most of them belonged to the PIPs subfamily. The signiﬁcantly upregulated LvAQP gene was Cluster-9981.114831 , and the signiﬁcantly downregulated LvAQP genes were Cluster-9981.112839 , Cluster-9981.107281 , and Cluster-9981.112777 . These genes might play a key role in responding to cold tolerance in the natural low-temperature growth stage of Ligustrum × vicaryi .


Introduction
Aquaporin is a protein located on the cytomembrane that controls the entry and exit of water in cells. Water uptake, transport across membranes and tissues are essential for plants growth and development, and the transmembrane transport of water molecules is mainly regulated by aquaporins. In biological membranes, plant aquaporins have a highly conserved Asn-Pro-Ala (NPA) motif structure, which plays a crucial role in the formation of water-selective channels [1]. It has been reported that AEFXXT motif located in the first helix (TM1) in plant aquaporins is highly conserved in almost all major intrinsic proteins (MIPs), but the exact function of the AEFXXT motif is still unclear [2]. The previous studies based on genomic data revealed that aquaporins constitute a huge gene family in plants. These aquaporins are divided into five main subfamilies according to their amino acid sequence [3]: plasma membrane intrinsic proteins (PIPs), tonoplast intrinsic proteins The seedlings overwintered naturally in the open field. The seedlings were sampled on the 24th of each month in September and November 2019 and in January and April 2020.

Determination of Cold Resistance of Ligustrum × vicaryi
The seedlings were taken at four sampling times. The seedlings were placed in an artificial climate chamber for low temperature treatment using the method of Di [38] ( Table 1). During the cooling period, the seedlings were kept at −3 • C for 5 h to keep the soil and air temperature consistent and were then continued to cool at the same temperature. Each set target temperature was maintained for 4 h. After that treatment, the seedlings were placed at 4 • C for 5 days and room temperature for 1 day, and the cold resistance of roots was measured by relative electrolyte leakage (REL).

. RNA Extraction and Detection
The seedlings were taken at four sampling times. The roots of the plants were washed by tap water to remove the soil, and the fine roots were rinsed with tap water, distilled water and ultrapure water in turn (the water was placed in Specimen Park in advance and the temperature was kept consistent with the environment), and then they were frozen in liquid nitrogen and stored in an ultra-low temperature freezer at −80 • C.
The RNA integrity was detected by agarose gel electrophoresis with 2% concentration, 150 V, 150 mA. The concentration of each RNA sample and its optical density in the wavelength range of 260 and 280 nm were measured by Nano Drop one spectrophotometer, and the OD 260 /OD 280 value was calculated to detect the purity of RNA; then the RNA was stored in an ultra-low temperature freezer at −80 • C.

cDNA Library Construction
First, magnetic beads with Oligo(dT) were used to enrich eukaryotic mRNA. Second, mRNA was broken into short fragments by adding fragmentation buffer. One-stranded cDNA was synthesized with six-base random hexamers using mRNA as a template. Third, double-stranded cDNA was formed by adding buffer, dNTPs, DNA polymerase I, and RNase H, which was purified by AMPure XP beads. The purified double-stranded cDNA was end-repaired and dA-tailed to ligate to sequencing connectors, and then fragment size selection was performed with AMPure XP beads. Finally, polymerase chain reaction (PCR) amplification was conducted, and the PCR products were purified with AMPure XP beads to obtain the final cDNA library. After the completion of the cDNA library construction, the initial quantification was operationalized by using Qubit 2.0, and then the library was diluted. Subsequently, the insert size of the library was tested. When the insert size met the expectation, the effective concentration of the library was accurately quantified by the Q-PCR method (effective library concentration > 2nM) to ensure the quality of the cDNA library.

RNA Data Analysis
The raw image data generated by the sequencer were transformed into raw data or raw reads by base calling. The results were stored in fastq format, which were part of the original file, including the sequence of reads and the sequencing quality of reads. Raw reads were processed to obtain clean reads by removing reads containing an adaptor, reads containing more than 10% of N and reads containing a small amount of low-quality sequences (the number of bases with quality value Q < 5 accounts for more than 50% of the entire reads).
The transcriptome data were assembled by Trinity v2.4.0 software with the following commands and parameters: Trinity-seq Type fq-max_memory 300G-left file_1. fq-right file_2.fq-CPU 50-full_clean up-KMER_SIZE 30-min_kmer_cov 5. Among the genes containing multiple transcripts, the sequence with the longest transcript was used as the basis for calculating expression, RSEM worked as the method for transcript abundance calculation, and trimmed mean of M-values (TMM) was used as the method for normalization between samples.
KEGG PATHWAY enrichment analysis on the results of variance analysis was performed by kobas software.

Identification of LvAQP Gene Family
Based on the transcriptome of Ligustrum × vicaryi and the study on the model plant Arabidopsis thaliana, the gene sequences of 35 Arabidopsis thaliana aquaporin genes (Table A1) were downloaded from NCBI (https://www.mhttps//www.ncbi.nlm.nih.gov/, accessed on 2 June 2021). The transcriptome database of Ligustrum × vicaryi was searched by blast homology retrieval method, and the LvAQP gene family was identified. The LvAQP genes were screened by MAFFT comparison software and a manual correction process using CD-HIT Suite (http://weizhong-lab.ucsd.edu/cdhit-web-server/cgi-bin/index.cgi?cmd= cd-hit-est, accessed on 10 June 2021) to remove redundancy. The screened LvAQP gene family protein sequences were analyzed for physicochemical properties such as number of amino acids, molecular weight, theoretical pI, aliphatic index and grand average of hydropathicity, using the online software ExPASy (https://web.expasy.org/protparam/, accessed on 25 June 2021).

Construction of Phylogenetic Tree
Multiple sequence alignment of candidate genes was performed by MAFFT software's E-INS-I strategy with necessary manual corrections. MEGA and PhyML3.0 software were used to construct the phylogenetic tree, and iTOL online software (https://itol.embl.de/ upload.cgi, accessed on 24 July 2021) was used for phylogenetic tree beautification. The constructed phylogenetic tree was analyzed by Alrt detection method and WAG model.

Conserved Motifs Analysis
LvAQP gene family was extracted using the MEME tool (https://meme-suite.org/ meme/tools/meme, accessed on 26 July 2021), with the following parameters: motif sequences, sites, width, E-value for each motif.

Gene Expression Profile Analysis
The LvAQP gene family expression profile was constructed using MeV software. The candidate genes were initially selected by differential gene expression analysis between natural low temperature treatment and non-low temperature treatment. Combining the constructed phylogenetic tree of LvAQP genes and the known cold resistance aquaporin genes of other plants aimed to find the homologs of known cold resistance genes of other plants within Ligustrum × vicaryi.  The UEIris II RT-PCR System for First-Strand cDNA Synthesis (with dsDNase) reverse transcription kit was used as follows: The RNA was denatured thermally at 65 • C for 5 min, immediately iced for more than 3 min, and then the reaction system was prepared as 20 µL of the reverse transcription system: 2 µL total RNA, 4 µL UEIris II RT MasterMix (5×), 13 µL RNase-free Water, 1 µL dsDNase.

LvAQP Gene
Reaction conditions: reverse transcription 37 • C for 2 min; 55 • C for 10 min; 85 • C for 10 s. After the reaction, it was stored at −20 • C.

Design of Primers for Quantitative Real-Time qRT-PCR
Nine LvAQP genes related to cold stress were selected for qRT-PCR. Primers were designed by Primer3 Plus and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China), and Ligustrum × vicaryi LvEF-1α was selected as an internal reference gene; the primer information is shown in Table S1.

qRT-PCR Reaction System and Reaction Conditions
According to the instructions of AugeGreen TM qPCR Master Mix reagent, Roche fluorescence quantitative PCR instrument lightcycl96 was used to detect the expression of target genes. Preparation of reaction solution for 20 µL reaction system: 10 µL 2× AugeGreen TM Master Mix, 7 µL ddH2O, 1 µL Forward Primer, 1 µL Reverse Primer, 1 µL cDNA template. qRT-PCR reaction procedure was as follows: 95 • C for 2 min; 40 cycles of 15 s at 95 • C and 60 s at 58 • C; 95 • C for 10 s; 65 • C for 60 s; 97 • C for 1 s.
The last step was to analyze the solubility curve of the amplification products to determine the specificity of the primers. There were 3 technical replicates and 3 biological replicates for each sample during qRT-PCR reactions. The results of qRT-PCR were calculated by the 2 −∆∆CT method to obtain gene expression.

Cold Resistance of Ligustrum × vicaryi during Natural Cold Stress Period
During the overwintering period, the cold resistance ability of Ligustrum × vicaryi gradually increased with the change of temperature in the early overwintering period, and gradually decreased in the late overwintering period ( Figure 1). The strongest cold resistance appeared in January, reaching −22.1 • C and 11.4 • C higher than that in September (−10.7 • C). The cold resistance in November was −15.2 • C, and the cold resistance in April was −13.0 • C, which was close to that in September.

LvAQP Gene Family Identification
Based on the Arabidopsis thaliana aquaporin gene family (Table A1), 58 candidate

LvAQP Gene Family Identification
Based on the Arabidopsis thaliana aquaporin gene family (Table A1), 58 candidate LvAQP genes were identified (Table 2). According to sequence alignment, the correlation of characteristic proteins and phylogenetic relationship, the 58 LvAQP genes were classified into four subfamilies: PIPs, TIPs, NIPs and SIPs, which contained 32, 11, 11 and 4 genes, respectively. The results of the analysis of the physicochemical properties of the LvAQP gene family proteins ( Table 2) showed that there were differences in the number of amino acids, molecular weight, theoretical pI, aliphatic index and grand average of hydropathicity of LvAQP gene family proteins. The number of amino acids in the protein sequence of the LvAQP gene family was 137-341. The molecular weight of the LvAQP gene family protein sequence was 13986.15-36029.94 Da. The theoretical pI of the LvAQP gene family protein sequences ranged from 5.34-9.91, and 40 of the 58 LvAQP proteins had a theoretical pI of greater than 7.5, indicating that most of the LvAQP gene family proteins are basic proteins. The aliphatic index of the LvAQP gene family protein sequences ranged from 89. 25-116.37. The grand average of hydropathicity of the LvAQP gene family proteins was positive, revealing that they were hydrophobic proteins.

Phylogenetic Analysis of LvAQP Gene Family
By constructing phylogenetic tree, the distribution and development of the 58 candidate aquaporin genes of the four subfamilies of the LvAQP gene family can be known ( Figure 2). The internal genes of the PIPs subfamily were more similar than that of the TIPs subfamily, NIPs subfamily, and SIPs subfamily. The PIPs subfamily had the largest genes, which were due to tandem repeats of some genes with similar structures on the chromosome. Among the four subfamilies, the PIPs subfamily contained seven pairs of tandem repeats genes, while only one tandem repeat gene in the NIPs and SIPs subfamilies ( Table 3). ( Figure 2). The internal genes of the PIPs subfamily were more similar than that of the TIPs subfamily, NIPs subfamily, and SIPs subfamily. The PIPs subfamily had the larges genes, which were due to tandem repeats of some genes with similar structures on the chromosome. Among the four subfamilies, the PIPs subfamily contained seven pairs o tandem repeats genes, while only one tandem repeat gene in the NIPs and SIPs subfami lies (Table 3).    Comparing the aquaporin gene family of Monocotyledonous Oryza sativa, Zea mays and Musa acuminata, and dicotyledonous Arabidopsis thaliana and Brassica rapa with that of Ligustrum × vicary, the results showed that the number of the PIPs gene subfamily was the largest, while the number of the SIPs gene subfamily was the smallest (Table 4). Furthermore, the distribution of the four gene subfamilies of the LvAQP gene family was generally the same as that of the subfamily members of other plants. Ligustrum × vicaryi is a dicotyledonous plant. Unlike other plants, in Ligustrum × vicaryi, the number of aquaporin genes in the PIPs subfamily was nearly two times higher than that of TIPs subfamily and NIPs subfamily, while it was similar to that of the TIPs and NIPs subfamilies in other plants. PIPs located on the plasma membrane were highly selective to the transport matrix, and they are critical for maintaining the water balance of cells in plants [39]. Thus, it isspeculated that the Ligustrum × vicaryi PIPs subfamily (LvPIPs) may play a major role in maintaining its own water balance of cells. In this study, a phylogenetic tree was constructed based on 35 Arabidopsis thaliana aquaporin genes, 35 Oryza sativa aquaporin genes, 58 candidate LvAQP genes and aquaporin genes related to cold stress in other plants ( Figure 3). Most of the aquaporin genes related to cold resistance were distributed in the PIPs subfamily (Figure 3), while the number of genes of the PIPs subfamily in Arabidopsis thaliana and Oryza sativa were relatively small (Table 4). There was only a pair of tandem repeats genes (At2G37170 and At2G37180) in the PIPs subfamily in Arabidopsis thaliana (Table A1), while there were seven pairs of tandem repeats genes in the LvPIPs subfamily. Therefore, the reason for the large number of LvPIPs may be that genes are relatively tightly distributed on chromosomes, and tandem duplication led to gene amplification.

LvAQP Gene Family Conserved Motifs Analysis
The identified LvAQP genes all contained conserved domains. Table 5 showes the LvAQP gene family containes 19 main conserved motifs. The distribution of 58 Lv conservative motifs isshown in Figure 4, and the four subfamilies share common servative motifs, such as motif1. Each subfamily has similar conserved sites, and the m bers of each subfamily contain similar conserved motifs, even the same, such as Clu 9981.115068 and Cluster-9981.109600 of the PIPs subfamily. Each subfamily contain own unique motifs. For example, all members of the PIPs subfamily contain motif 3, m 4, motif 5, and motif 11. All members of the TIPs subfamily contain motif 7 and mot All members of the NIPs subfamily contain motif 12, motif 14, and motif 15. The conn motif 17 and motif 1 are in the SIPs subfamily. Each LvAQP subfamily was highly served during the process of evolution, which was beneficial to the phylogeny o LvAQP gene family.

LvAQP Gene Family Conserved Motifs Analysis
The identified LvAQP genes all contained conserved domains. Table 5 showes that the LvAQP gene family containes 19 main conserved motifs. The distribution of 58 LvAQP conservative motifs isshown in Figure 4, and the four subfamilies share common conservative motifs, such as motif1. Each subfamily has similar conserved sites, and the members of each subfamily contain similar conserved motifs, even the same, such as Cluster-9981.115068 and Cluster-9981.109600 of the PIPs subfamily. Each subfamily contains its own unique motifs. For example, all members of the PIPs subfamily contain motif 3, motif 4, motif 5, and motif 11. All members of the TIPs subfamily contain motif 7 and motif 19. All members of the NIPs subfamily contain motif 12, motif 14, and motif 15. The connected motif 17 and motif 1 are in the SIPs subfamily. Each LvAQP subfamily was highly conserved during the process of evolution, which was beneficial to the phylogeny of the LvAQP gene family.   All the known aquaporin genes related to cold stress contain several common gene sequence fragments ( Figure 5), namely IAEFXXT, GIAW, GGMI, LVYCTAG, SGGHINPAVT, GTFVLVYTVF and ATD, which may play a key role in resisting cold. The above fragments in LvAQP came from motif 2, motif 3, motif 4, and motif 6. Most of these motifs were dis-tributed in the PIPs subfamily of LvAQP. Therefore, the PIPs subfamily might be important for Ligustrum × vicaryi under cold stress.
All the known aquaporin genes related to cold stress contain several common gene sequence fragments (Figure 5), namely IAEFXXT, GIAW, GGMI, LVYCTAG, SGGHINPAVT, GTFVLVYTVF and ATD, which may play a key role in resisting cold. The above fragments in LvAQP came from motif 2, motif 3, motif 4, and motif 6. Most of these motifs were distributed in the PIPs subfamily of LvAQP. Therefore, the PIPs subfamily might be important for Ligustrum × vicaryi under cold stress.

Analysis of LvAQP Gene Expression Pattern
The transcript abundance of LvAQP was analyzed in four sampling times, and combining the phylogenetic relationship between cold stress aquaporin genes of various plants and of LvAQP genes (Figure 3) was helpful in identifying the specific expression patterns of individual genes of the LvAQP gene family.
The expression of 58 LvAQP genes changed in September, November, January and April ( Figure 6). Transcriptional analysis showed that the PIPs subfamily and TIPs subfamily contained a relatively high expression in four sampling times. Compared to September, 8% of LvAQP gene expression increased in November and January and decreased in April; 21% of LvAQP gene expression decreased in November and January and increased in April; 24% of LvAQP gene expression increased in November, decreased in January, and increased in April; 17% of LvAQP gene expression decreased in November, increased in January, and decreased in April; 21% of LvAQP gene expression increased in November and decreased in January and April; 3.4% of LvAQP gene expression decreased in November and increased in January and April; and 5% of LvAQP gene expression decreased consecutively in November, January, and April. According to relevant research, the overexpression of MusaPIP1;2 in Musa acuminata enhanced plant freezing resistance [40]; in Arabidopsis thaliana, the overexpression of AtPIP1;4 and AtPIP2;5, along with repressed expression of other PIPs family members enhanced plant cold resistance [25]; in Oryza sativa, there was increased expression of OsPIP2;5 and OsPIP2;7 and decreased expression of OsPIP1;3, which helped to improve cold resistance [28,41]. Studies had shown that plants enhance cold resistance by overexpressing or inhibiting the expression of aquaporin genes under cold stress. Therefore, in this study, the researchers selected LvAQP genes, whose gene expression increased in November and January and decreased in April, and whose gene expression decreased in November and January and increased in April in four sampling times, as the target genes.

Analysis of LvAQP Gene Expression Pattern
The transcript abundance of LvAQP was analyzed in four sampling times, and combining the phylogenetic relationship between cold stress aquaporin genes of various plants and of LvAQP genes (Figure 3) was helpful in identifying the specific expression patterns of individual genes of the LvAQP gene family.
The expression of 58 LvAQP genes changed in September, November, January and April ( Figure 6). Transcriptional analysis showed that the PIPs subfamily and TIPs subfamily contained a relatively high expression in four sampling times. Compared to September, 8% of LvAQP gene expression increased in November and January and decreased in April; 21% of LvAQP gene expression decreased in November and January and increased in April; 24% of LvAQP gene expression increased in November, decreased in January, and increased in April; 17% of LvAQP gene expression decreased in November, increased in January, and decreased in April; 21% of LvAQP gene expression increased in November and decreased in January and April; 3.4% of LvAQP gene expression decreased in November and increased in January and April; and 5% of LvAQP gene expression decreased consecutively in November, January, and April. According to relevant research, the overexpression of MusaPIP1;2 in Musa acuminata enhanced plant freezing resistance [40]; in Arabidopsis thaliana, the overexpression of AtPIP1;4 and AtPIP2;5, along with repressed expression of other PIPs family members enhanced plant cold resistance [25]; in Oryza sativa, there was increased expression of OsPIP2;5 and OsPIP2;7 and decreased expression of OsPIP1;3, which helped to improve cold resistance [28,41]. Studies had shown that plants enhance cold resistance by overexpressing or inhibiting the expression of aquaporin genes under cold stress. Therefore, in this study, the researchers selected LvAQP genes, whose gene expression increased in November and January and decreased in April, and whose gene expression decreased in November and January and increased in April in four sampling times, as the target genes. By analyzing the relative transcript abundance profile of LvAQP genes ( Figure 6) and the phylogenetic relationship between cold stress aquaporin genes of various plants and LvAQP genes (Figure 3), 20 LvAQP genes that responded to cold stress were determined: Cluster-9981.109600, Cluster-9981.112839, Cluster-9981.112265, Cluster-9981 8803. Among the determined 20 LvAQP genes, 11 genes were part of the PIPs subfamily, fivegenes were part of the TIPs subfamily, and two genes were part of the NIPs subfamily and SIPs subfamily, separately.
Among the 20 LvAQP genes identified in response to cold stress, the expression of Cluster-9981.114831 was significantly up-regulated during the two periods of lowest natural temperature in November and January and the most cold-resistant period in January, while the expression of three genes was significantly downregulated, namely, Cluster-9981.112839, Cluster-9981.107281 and Cluster-9981.112777. All the significantly upregulated genes contained motif 6, and all the significantly downregulated genes contained motif 1 and motif 2, which were basically consistent with the common special motifs reported in aquaporin genes related to cold stress. It was speculated that the key role of some AQP genes in Ligustrum × vicaryi for cold resistance might respond to the presence of these specific modular motifs. Among the 20 LvAQP genes identified in response to cold stress, the expression of Cluster-9981.114831 was significantly up-regulated during the two periods of lowest natural temperature in November and January and the most cold-resistant period in January, while the expression of three genes was significantly downregulated, namely, Cluster-9981.112839, Cluster-9981.107281 and Cluster-9981.112777. All the significantly upregulated genes contained motif 6, and all the significantly downregulated genes contained motif 1 and motif 2, which were basically consistent with the common special motifs reported in aquaporin genes related to cold stress. It was speculated that the key role of some AQP genes in Ligustrum × vicaryi for cold resistance might respond to the presence of these specific modular motifs.

KEGG Enrichment Analysis of Differentially Expressed Genes
KEGG pathway enrichment analysis of differentially expressed genes was conducted under natural cold stress in Ligustrum × vicaryi (Table 6). A total of 12,872 differentially expressed genes were distributed in 338 pathways, and 10 of them showed significant differences (p < 0.05). The differentially expressed genes were significantly enriched in ribosome (ko03010), starch and sucrose metabolism (ko00500), plant hormone signal transduction (ko04075). Antigen processing and presentation 2.38E-02 72 ko00520 Amino sugar and nucleotide sugar metabolism 3.49E-02 122

Expression Verification of Cold-Responsive LvAQP Target Ggenes
The nine screened LvAQP genes that target differentially expressed gene were verified by real-time PCR. The qRT-PCR results were first calculated by the 2 − CT method, followed by log calculation based 2. The change of log 2 multiples for the real-time fluorescence quantification of the nine target genes areshown in Table 7. Although the real-time PCR results of individual genes deviated from the RNA-Seq results in terms of differential fold, the up-regulated and down-regulated expression trends between them were consistent (Figure 7). In addition, the correlation analysis between the results of the qRT-PCR analysis and RNA-seq sequencing results showed that the correlation coefficient R 2 reached 0.70 (Figure 8), indicating that the transcriptome sequencing results of Ligustrum × vicaryi cold resistance were reliable.

Discussion
The number of genes encoding aquaporin of Ligustrum × vicaryi was more than that in Arabidopsis thaliana, especially in the PIPs subfamily due to gene amplification. In this study, 58 candidate LvAQP genes were found. Phylogenetic analysis showed that these 58 LvAQP genes can be divided into four subfamilies: PIPs, TIPs, NIPs and SIPs. A fifth subfamily has also been reported: XIPs, which is a class of atypical non-specific intrinsic

Discussion
The number of genes encoding aquaporin of Ligustrum × vicaryi was more than that in Arabidopsis thaliana, especially in the PIPs subfamily due to gene amplification. In this study, 58 candidate LvAQP genes were found. Phylogenetic analysis showed that these 58 LvAQP genes can be divided into four subfamilies: PIPs, TIPs, NIPs and SIPs. A fifth subfamily has also been reported: XIPs, which is a class of atypical non-specific intrinsic

Discussion
The number of genes encoding aquaporin of Ligustrum × vicaryi was more than that in Arabidopsis thaliana, especially in the PIPs subfamily due to gene amplification. In this study, 58 candidate LvAQP genes were found. Phylogenetic analysis showed that these 58 LvAQP genes can be divided into four subfamilies: PIPs, TIPs, NIPs and SIPs. A fifth subfamily has also been reported: XIPs, which is a class of atypical non-specific intrinsic aquaporins. It was absent in Arabidopsis thaliana, Oryza sativa, Zea mays and Ligustrum× vicaryi. Plasma membrane intrinsic proteins are highly selective for the transporting matrix, and they play an important role in maintaining cell water balance under various adversities [15]. Studies have shown that plants resist abiotic stress by regulating the expression and activity of PIPs in the plasma membrane [39,[42][43][44]. Plants mainly regulate their response to stress through the expression or inhibition of PIPs genes of the aquaporin family. Under natural low temperature adversity, maintaining water balance in the body is a considerable challenge to Ligustrum × vicaryi. At this time, the transmembrane transport of water in Ligustrum × vicaryi mainly depends on the PIPs subfamily of the aquaporin family. This study found that the number of PIPs subfamily was the largest in the LvAQP gene family, which was consistent with the results of previous studies [4][5][6][7][8][9][10]. The differences from previous studies were that Arabidopsis thaliana and Oryza sativa have 13 and 11 PIPs genes, respectively, while this study found 32 PIPs subfamily genes in Ligustrum × vicaryi. The number of genes in the LvPIPs subfamily was much higher than that of other plants. When a certain gene family has obvious gene clusters on the chromosome, it is often accompanied by the gene expansion mechanism of tandem replication [45]. The large number of the PIPs subfamilies of LvAQP gene family was caused by the expansion and tandem duplication of some genes with similar structure in the gene cluster. In this study, the number of LvPIPs genes was higher than that of other plants. There were 11 of the 20 aquaporin genes screened that were related to low temperature stress belonged to the PIPs subfamily. The result of 11 genes belonged to the PIPs subfamily was in accordance with previous studies on aquaporins in response to cold stress, suggesting that the PIPs subfamily of aquaporin might play a major role in responses to cold stress in Ligustrum × vicaryi [11,[28][29][30]32,33]. Unlike in previous studies, two genes of the SIPs subfamily in the LvAQP gene family also responded to cold stress.
In the face of cold stress, plants generally respond to stress by regulating water homeostasis in the body, in which aquaporin proteins are one of the key pathways of water transport [46][47][48][49]. The expression patterns of aquaporins in various plant tissues are different, which indicates that aquaporins may have different functions in plants [50]. After freezing treatment, the low-temperature-tolerant Zea mays variety z7 maintained root hydraulic conductivity and water transport by expressing a large amount of aquaporins to reduce freezing damage [51]. The aquaporins PIP1 and PIP2 of Arabidopsis thaliana cooperated synergistically in the roots under cold stress to affect root hydraulic conductivity and to regulate plant cold resistance [52]. Overexpression of PtPIP2;5, PtPIP2;1 and PtPIP2;3 in Populus trichocarpa affected its response to cold stress and osmotic stress [53]. Under cold stress, the overexpression of banana MaPIP2;7 lowered the MDA content and electrolyte leakage in the plant, while the content of chlorophyll, proline, soluble sugar and ABA was higher, thereby enhancing the tolerance to various stresses such as the cold [54]. The overexpression of MaSIP2;1, OsPIP2;7, and TaAQP7 (PIP2) regulated the osmotic balance in plants, reduced membrane damage and oxidation, and adjusted the levels of hormones such as ABA and GA to improve the cold tolerance of plants [30,36,43].
In this study, the phylogenetic comparison between LvAQP genes and reported aquaporin genes related to cold stress in other plants as well as the changes of aquaporin genes transcription abundance in four sampling times was conducted to identify the specific expression patterns of individual genes of the gene family under natural cold stress. The 20 aquaporin genes that responded to cold stress were screened from the 58 LvAQP genes; eleven belonged to the LvPIPs subfamily, five belonged to the LvTIPs subfamily, and two belonged to the LvNIPs subfamily and LvSIPs subfamily, separately, which indicated that genes of the PIPs subfamily played a major role in response to natural cold stress. Among these 20 aquaporin genes of Ligustrum × vicaryi that responded to cold stress, all the significantly upregulated genes contained motif 6, while all the significantly downregulated genes contained motif 1 and motif 2. It was speculated that motif 1, motif 2 and motif 6 might play an important role in response to cold stress when Ligustrum × vicaryi is under a natural low temperature. In the analysis of the reported cold stress-related AQP gene sequences of other plants, we found that the gene sequence SGGHINPAVT was present in motif 2 and GIAW and GGMI were present in motif 6; thus, it was further speculated that the gene sequences of SGGHINPAVT, GIAW and GGMI might play a major role in the response to cold stress in Ligustrum × vicaryi. However, the AEFXXT motif, which was conserved in almost all MIPs in previous studies, was not conserved in all significantly upregulated and significantly downregulated genes in Ligustrum × vicaryi in response to cold stress. Therefore, we speculated that the AEFXXT motif might not be the key motif in genes responding to cold stress. From the determination of cold resistance of Ligustrum × vicaryi, it can be seen that Ligustrum × vicaryi was most resistant to cold in January during the natural overwintering process, and the cold resistance of the plant changed with the change of time. In this study, 75% of the LvAQP genes that were significantly related to cold stress decreased in November and January, and their expression increased in April, which is consistent with the results of transcriptome analysis of Arabidopsis thaliana, Oryza sativa, and the roots and leaves of Zea mays [25,43,51]. In winter, low temperatures can easily lead to freeze thaw embolism of plants, which blocks water transport and leads to withering. At this time, aquaporin may be involved in embolization repair. Low soil temperature limits the absorption of water by roots, leading to a water imbalance. Low soil temperature can reduce or increase the activity of aquaporin in roots, but appropriate low temperature acclimation can promote the abundance of AQP in roots. In the process of natural cold stress, with the enhancement of cold resistance, Ligustrum × vicaryi regulated the decrease or increase in the expression of aquaporin genes and the corresponding protein activity, and adjusted root hydraulic conductivity, thus maintaining the water balance in the plant, resisting the effects of natural low-temperature stress, and ensuring normal life activities.
Aquaporins are important membrane functional proteins in many physiological reactions, which play a key role mainly through transcriptional regulation, post-translational modification and subcellular localization [55,56]. Plasma membrane intrinsic proteins and tonoplast intrinsic proteins are located on the inner chloroplast membrane and thylakoid membrane [3]. KEGG enrichment analysis of Ligustrum × vicaryi genes showed that they responded to cold stress mainly through the sucrose metabolism pathway and plant hormone signal transduction pathway. It was speculated that some genes of the PIPs and TIPs subfamilies on the plasma membrane and in the chloroplast were upregulated or downregulated, which would enhance the cold resistance of Ligustrum × vicaryi by regulating the synthesis and transformation of soluble sugar or starch. After feeling a natural low temperature, differentially expressed genes related to hormone signaling were enriched, and pathways such as ABA signaling were turned on under low temperature stress, thus inducing the expression of downstream regulatory genes. Then, the expression of AQP genes changed in order to regulate the synthesis of corresponding proteins and other macromolecules, to stabilize the membrane structure, and to reduce the water transport rate to avoid low temperature damage of Ligustrum × vicaryi.

Conclusions
In this research, the gene expression of LvAQP under natural cold stress was studied. We identified 58 candidate LvAQP genes. Based on phylogenetic analysis, the 58 candidate LvAQP genes were divided into four subfamilies: 32 belonged to the PIPs subfamily, 11 belonged to the TIPs subfamily, 11 belonged to the NIPs subfamily and 4 belonged to the SIPs subfamily. The number of genes in the PIPs subfamily was nearly twice as large as that in other plants. The LvAQP gene family contained nine pairs of tandem repeats genes, which had high conservatism in the process of evolution by searching for conserved motifs. We obtained 20 differentially expressed LvAQP genes under natural cold stress. Among the 20 differentially expressed genes, 11 belonged to the LvPIPs subfamily. Among the 20 differentially expressed genes, the significantly up-regulated gene was Cluster-9981.114831; the significantly down-regulated genes were Cluster-9981.112839, Cluster-9981.107281 and Cluster-9981.112777. These four LvAQP genes might play important roles in response to low temperature stress. The results laid the foundation for further exploration of cold resistant aquaporin genes and biological function verification of Ligustrum × vicaryi.