Genome-wide analysis of the aquaporin genes in melon (Cucumis melo L.)

Melon (Cucumis melo L.) is a very important crop throughout the world and has great economic importance, in part due to its nutritional properties. It prefers well-drained soil with low acidity and has a strong demand for water during fruit set. Therefore, a correct water balance—involving aquaporins—is necessary to maintain the plants in optimal condition. This manuscript describes the identification and comparative analysis of the complete set of aquaporins in melon. 31 aquaporin genes were identified, classified and analysed according to the evolutionary relationship of melon with related plant species. The individual role of each aquaporin in the transport of water, ions and small molecules was discussed. Finally, qPCR revealed that almost all melon aquaporins in roots and leaves were constitutively expressed. However, the high variations in expression among them point to different roles in water and solute transport, providing important features as that CmPIP1;1 is the predominant isoform and CmTIP1;1 is revealed as the most important osmoregulator in the tonoplast under optimal conditions. The results of this work pointing to the physiological importance of each individual aquaporin of melon opening a field of knowledge that deserves to be investigated.

Melon (Cucumis melo L.) is a eudicot diploid species (2n = 24) and one of the most important cucurbits. The great importance of melon is due not only to its great economic importance but also to its nutritional properties as the presence of health-promoters including β-carotenes, folic acid, phenolic acids, vitamins A and C and minerals 16,17 , being of special interest in nutrition and human healthcare 18 . The crop highly demand an appropriate water balance 19 related with its capacity to adapt quickly under abiotic stress conditions. This suggests strong control of water and nutrients transport in membranes, pointing to aquaporins as one of the most interesting targets in these plants.
Therefore, as each aquaporin is structurally unique and that simple variations of specific residues can produce an altered solute selectivity, the aim of this study was to identify, compare and analyse the complete set of aquaporins in melon. The comparative analysis of aquaporins was performed to understand the evolutionary relationship of melon with related plant species, to predict their possible roles in the transport of other compounds apart from water and to determine their constitutive expression (determining structural importance of NPA motifs, the ar/R region and FPs). Analysis by qPCR was carried out in roots and leaves after the design of primers and a comparison with RNA-seq results was done with existing databases.

Material and methods
Identification of putative C. melo aquaporins (CmAQPs). The complete set of aquaporin protein sequences of watermelon (Citrullus lanatus L.) 20 , identified by Zhou et al., 2019, was used as a template in the PSI-Blast tool (https ://blast .ncbi.nlm.nih.gov/Blast .cgi) with default parameters and the three-iterations method. This was performed against the C. melo database and using the non-redundant protein sequences (nr) database.
For the phylogenetic analysis and tree construction, C. melo aquaporins protein sequences were aligned with all the characterised aquaporins from C. lanatus 20 , A. thaliana 2 and Cucumis sativus L. 47 (Fig. 1). Mega X was employed to build the phylogenetic tree. We aligned all the sequences with the MUSCLE method and then used a Neighbour Joining (NJ) algorithm employing 1000 bootstrap replicates, a Poisson model and pairwise deletion. The same analysis and tree construction were performed for the C. melo, A. thaliana, Z. mays 48 and O. sativa 49 phylogenetic tree (Supplementary Figure S2). Plant material and growth conditions. Seeds of melon, C. melo var. Grand Riado, supplied by SAKATA SEED IBERICA S.L.U, were hydrated with deionised water and aerated for 16 h. After this, they were germinated in vermiculite, in the dark at 28 °C, for 2 d. They were then transferred to a controlled-environment chamber having a 16-h light and 8-h dark cycle with temperatures of 25 and 20 °C and relative humidity of 60 and 80%, respectively. Photosynthetically active radiation (PAR) of 400 µmol m −2 s −1 was provided by a combination of fluorescent tubes (PHILIPS TLD 36 W/83, Jena, Germany and Sylvania F36 W/GRO, Manchester, NH, USA) and metal halide lamps (OSRAM HQI, T 400 W, Berlin, Germany). After 3 d, once the seeds germinated, 16 seedlings were placed in 15-L containers, 4 seedlings per container, with continuously-aerated Hoagland nutrient solution 50 . The solution was replaced completely every week. After a month of growth, 6 plants were harvested separating roots and leaves from the rest of the plant, instantly frozen in liquid nitrogen and stored at − 80 ºC. The experimental design was completed randomized design (CRD).
Primers design. The primer sets used to amplify each aquaporin gene were specifically designed in the 3′ or 5′ non-coding region of each gene, in order to avoid the non-specific amplification of other aquaporin genes, since in the coding region the homology between their sequences is very high 51 . For this same reason, the primer design was carried out manually, meeting the specific requirements imposed by the high homology of the sequences used and the technique used for the analysis. Virtual analysis of melting temperature, primer hairpins, self-dimers, hetero-dimers and individual and total ΔG was performed with PCR Primer Stats (https ://www. bioin forma tics.org/sms2/pcr_prime r_st) and IDT Oligo Analyzer Tools (https ://eu.idtdn a.com/calc/analy zer/). The ΔG accepted for dimer analysis was always less than − 6.5 kcal/mol ( Table 1). The specificity of the amplicons was checked using the virtual nucleotide basic local alignment search tool (NCBI nucleotide BLAST: https :// blast .ncbi.nlm.nih.gov/Blast .cgi) and then physically, with the total DNA extracted from 50 mg of frozen sample C. lanatus (stars) and C. sativus (squares). We used MUSLE to align the protein sequences and the NJ method (with 1000 bootstrap replications) to build the tree, all with MEGA X. Phylogenetic tree design has been done with the online tool "Interactive Tree Of Life" (iTOL; https ://itol.embl.de/). The different MIP sub-families are highlighted as: PIPs (green; PIP1s light green and PIP2s dark green), TIPs (sky blue), NIPs (orange), SIPs (pink) and XIPs (purple). The abbreviation of the species is as follows: Cm (C. melo), At (A. thaliana), Cl (C. lanatus), Cs (C. sativus). The meanings of the suffixes in the C. sativus names are: lk: like, pb: probable, pd: predicted.  F:TGC ACC CTT TGG ACT  TCT TC  R:GAT GTA GTA GCC ATC  CCG TAAAC   104  ---98.7  ---CmPIP1;1   F: GCT TCT TCA ATC AAT  CAC AGC  R: ACC ATT ACA TAG CTT  CAT 53 , were checked in each cDNA used in the quantitative PCR quantification (qPCR) and were measured using a Visual basic application for Excel (GeNorm) that automatically calculates the gene stability 54 . CmRAN (encoding the GTP-binding nuclear protein) was then selected as the reference gene for the standardisation of each sample. The sequences and features of the primers used for the 31 melon aquaporin genes and one constitutively-expressed gene are shown in Table 1.

Quantitative real-time PCR (RT-qPCR) analyses.
To compare the expression of all aquaporin genes, RT-qPCR was carried out. It was performed using 2 μL of 1:2, 1:5 or 1:10 diluted cDNA samples, depending on the gene being analyzed (see Table 1), in 10 μL of reaction medium containing 500 nM gene-specific primers and www.nature.com/scientificreports/ 5 μL of Power SYBR Green PCR Master Mix (Applied Biosystems by Thermo Fisher Scientific), in a QuantStudio 5 Flex, a Real-Time qPCR system (Applied Biosystems by Thermo Fisher Scientific) following the manufacturer's instructions. The qPCR program consisted of a 10 min initial denaturation at 95 °C and then amplification in a two-step procedure: 15 s of denaturation at 95 °C and 60 s of annealing and extension at a primer-specific temperature for 40 cycles, followed by a dissociation stage. Data collection was carried out at the end of each round in step 2. These conditions were used for both target and reference genes, and the absence of primer-dimers was checked in controls lacking templates. Real-time PCR analysis was performed on 3-6 independent samples for each treatment (biological replicates) and each sample reaction was carried out in triplicate (technical replicates) in 96-well plates. The transcript levels were calculated using the 2−ΔΔCt method 55 . Negative controls without cDNA were used in all PCR reactions. Finally, the normalised expression levels were rescaled and presented as relative units (ru) with respect to the most highly-expressed aquaporin, CmTIP1;1, which was assigned a value of 100. The extremely high expression of CmTIP1;1 masks the presence of other aquaporins, which is why we separated them into three different groups depending on their relative levels of expression.
Comparison with previous RNA-seq. Our qPCR analyses were compared with the RNA-seq analyses in databases, namely, the work of Latrasse et al. (2017) 56 . For this, the levels of all RT-qPCRs were relativized to that of the gene that had the highest expression levels, as previously described. The same procedure was followed with all the genes analysed by RNA-seq. In this case, the maximum value was for PIP1;1 in roots, which received the same fixed value of 100, and the rest were relativized to it. Finally, the values of each aquaporin in the RNAseq were compared with those obtained from our RT-qPCR. After being relativized, all aquaporins data were separated in groups of similar relative expression, as in the previous case of the qPCR data. Ethical approval. This article does not contain any studies with human or animal subjects.

Consent to participate.
All the authors referred in this study acknowledge the authorship and agreed with the content.

Consent of publication.
All the authors referred in this study give explicit consent to submit this article and consent to the publication of the data presented here.

Results
Genome-wide identification of CmAQP genes. A search of the whole genome for aquaporins proteins revealed 57 matches in the NCBI protein database, 34 corresponding to C. melo, while the rest of the nonselected sequences were specific to concrete varieties. First, we analysed each matched sequence and searched for its accession in the Cucurbit genomics database. During this procedure we found three sequences (MELO3C025166.2, MELO3C000776.2 and MELO3C025165.2) that were not previously located in our first search of the NCBI protein database. After this, we analysed all the sequences, including these new finds, which resulted in the elimination of some sequences based on the following. The CmTIP1;3 gene had three transcript variants in the nucleotide database and only one of them had the characteristic NPA motifs when transcribed to protein. So, the other two transcript variants were omitted.
Also, we found a high nucleotide sequence similarity among five sequences, the three new sequences found in the Cucurbit genomics database and two sequences from the 34 found in the NCBI database (XP_008466925, XP_008463222.1). Only XP_008463222.1 (from now on named CmPIP1;1) was a complete protein sequence including the two characteristic NPA motifs. Sequence alignments showed that XP_008466925 was, indeed, a partial sequence of CmPIP1;1 with 100% sequence identity. In addition, the other three sequences also had a high similarity to CmPIP1;1 but the alignment seemed to be shifted, so we thought these sequences were from the upstream region of the CmPIP1;1 gene, integrating only the second NPA motif. To prove this theory, we extracted the sequence downstream and upstream of the CmPIP1;1 gene in chromosome 10 and then aligned it against these three unknown sequences. As a result, all the sequences had at least 98-99% identity with the sequence including the CmPIP1;1 gene and the alignments were located upstream of CmPIP1;1, containing only the second NPA motif in their sequences, as we thought at the beginning. With these results and due to the lack of the two NPA motifs, we decided to omit the sequences from further analysis. The presence of this type of sequence in the database could be due to a non-curated identification method without experimental evidence.
In all the analyses carried out, a total of 31 putative aquaporin genes were identified (Table 2).
Nomenclature and classification. For the 31 putative aquaporin genes finally selected, their amino acid sequences were aligned for a phylogenetic analysis that divided them into sub-families, to help us to name them correctly. The C. melo aquaporins were named relative to their homology and phylogenetic relationships with those of C. sativus, C. lanatus and A. thaliana (Fig. 1). The proteins were grouped into five sub-families consisting of twelve PIPs, nine TIPs, eight NIPs, two SIPs and one XIP (Fig. 1). The PIP sub-family was divided into two groups (PIP1 and PIP2), with two members in the www.nature.com/scientificreports/ PIP1 group and ten in PIP2. The TIP sub-family was divided into five groups: TIP1, with three members, TIP2, with two members and one member each in groups TIP3, TIP4 and TIP5. The NIPs were divided into NIP1, NIP4 and NIP7, with one member each, and groups NIP2 and NIP5, with two members each. The SIP sub-family was formed by SIP1 and SIP2, with one member each. Lastly, only one sequence represented the XIP group.
Chromosomal location, protein features and subcellular localisation prediction. These 31 putative aquaporin-encoding genes were located on all chromosomes except chromosome 2, in a non-uniform manner: chromosome 9 incorporated six aquaporins genes, while chromosomes 4 and 5 had five genes each, chromosome 8 had three genes, two genes were identified on each of chromosomes 1, 3, 6, 10 and 11 and chromosomes 7 and 12 possessed only one putative aquaporin gene (Fig. 2). Exon-intron analysis showed a similar distribution on exons structures among each aquaporin family displaying 3-5 exons per gene (Supplementary Figure S1). The list of identified aquaporins, together with basic statistics on the primary protein sequences, is reported in Table 2. The PIPs protein length ranged from 276 (CmPIP2;5) to 292 residues (CmPIP1;1), with a median of 283 amino acids (aa) per protein, while that of the TIP sub-family ranged from 247 to 284 residues (CmTIP4;1 and CmTIP3;1, respectively), with a median of 255 aa. In the case of the NIPs, the length varied between 250 (CmNIP5;2) and 305 residues (CmNIP6;1), with an average of 276 aa. Finally, the SIP sub-family (with only two members) had a mean length of 275 aa. Motifs structure of each AQP of the same subfamily showed similar motifs structures (PIPs, TIPs and NIPs), NPA motif, and a large motif (YRALIAEFIATLLFLFVGVLT-VIGYSKQTDTASLGGIG) conserved in all the sequences. However, SIPs and XIP1.1 showed different motif distribution, with some similarities between themselves. (Supplementary Figure S3). Table 2 also displays the Pi and Mw data. The Pi varied among the families but was more or less constant among the sub-family members, with some exceptions. The PIPs sub-family had an average Pi of 8.58. In contrast, the TIPs group showed a mean Pi of 6.08, the lowest of all the sub-families, but CmTIP5.1 stood out with a Pi of 8.31, quite high for this group. The median Pi was 7.81 in the NIPs and, lastly, the SIPs displayed the highest Pi of all the sub-families, with an average of 9.63. The Mw did not vary much among the sub-families, unlike the Pi. The PIPs, NIPs and SIPs had mean values of 30.27 kDa, 29.09 kDa and 29.32 kDa, respectively; only the TIPs presented a notable difference, with an average of 26.39 kDa.
One of the principal features that define not only the aquaporins but also Integral Membrane Proteins (IMPs) is the presence of transmembrane helices. In the case of aquaporins, they must number at least six. We present, in Table 2, a prediction of the possible number of transmembrane motifs of each protein. For four of them, only five transmembrane helices were predicted and in order to confirm these results, we analysed the 3D structures of these proteins (Fig. 3); it can be seen that these four proteins displayed six transmembrane helices.
Prediction of the subcellular localisation, based on bioinformatics tools (the Plant-mPLoc and Wolf-PSORT programs), predicted all the PIP, SIP and XIP aquaporins situated on the plasma membrane. Most of the TIPs     Table 3. The AQP phylogenetic framework can be used to predict the putative function of individual AQPs on the basis of orthologous genes from A. thaliana and other species 11 . Thus, the identification of conserved motifs in each subfamily and in each cluster of orthologous genes offers a framework for studying their possible functional implication 57 . Based on this principle, a comparative analysis of the main important residues of orthologous genes in C. melo and A. thaliana and their homologues in other species such as rice and maize was performed (Supplementary Table S1), as the transport of different solutes by their aquaporins has been studied previously, both functionally and experimentally. In addition to this comparison of homologues, mutational studies, pore structure analyses and experimental studies were considered, to propose the transport possibilities of each C. melo aquaporin, specified in Table 3. www.nature.com/scientificreports/ Expression analysis. After primer verification, the expression of 31 aquaporins of C. melo was analysed in root and leaf tissue by RT-qPCR. All aquaporins were detected in both tissues, except PIP2;9, which was only detected in roots, and XIP1;1, that could not be detected in any sample. Figure 4 shows the aquaporins grouped according to their levels of expression (see Material and Methods). Group 1 contains the aquaporins that showed higher expression (10-100 ru) in both tissues: TIP1;1, PIP1;2 and PIP1;1 (Fig. 4a). Group 2 contains those with medium levels of expression (1-10 ru): PIP2;2, PIP2;3, PIP2;6, PIP2;10, TIP3;1, NIP2;1, NIP2;2, NIP5;1 and NIP5;2 (Fig. 4b). Finally, group 3 contains the aquaporins PIP2;1, PIP2;4, PIP2;5, PIP2;7, PIP2;8, PIP2;9, TIP1;2, TIP1;3, TIP2;1, TIP2;2, TIP4;1, TIP5;1 NIP1;1, NIP4;1, NIP6;1, NIP7;1, SIP1;1, SIP1;2 and XIP1;1, which showed low levels of expression (0-1 ru) (Fig. 4c).
The aquaporin that showed the highest gene expression in the analysed tissues was TIP1;1. Among the PIP subgroups, the PIP1s clearly showed higher expression levels, PIP1;1 being stronger in the leaves and both PIP1;1 and PIP1;2 in the roots, while of the PIP2s the most-highly expressed was PIP2;10, in both roots and leaves, this latter tissue also showing high expression of PIP2;6. Nor should the presence of PIP2;2 and PIP2;3, in both tissues, and of PIP2;4, in roots, be neglected. Within the TIPs subfamily, TIP1;1 exhibited the highest expression in both roots and (especially) leaves, followed by TIP3;1, in both tissues. Attending to the NIPs groups, the most expressed were NIP2;1, NIP2;2, NIP5;1 and NIP5;2; the expression of NIP2;1 was more important in leaves while both NIPs5 and NIP2;2 were expressed in roots and leaves at similar levels. Regarding the rest of the aquaporins, it should be noted that most of them appeared to a greater extent in the roots, being either lower, or practically nil, their presence in leaves, with no significant differences between them. As exceptions, PIP2;5, TIP2;2 and TIP1;3 had practically the same level of expression in both tissues, although slightly higher in leaves.

Discussion
The availability of the whole genome sequence of melon would facilitate genome-wide analysis to identify the complete set of aquaporins. Also, some analyses of the expression of diverse aquaporins with RNA-seq techniques have been performed; thanks to these it is possible to find in the various databases the nucleotide sequences of a large number of melon aquaporins genes, although many of them come from diverse origins and have not been ordered, compared and classified. In this way, the number of aquaporins genes identified in melon (31) was comparable with the number found in other plant species, such as rice (33) 49 , A. thaliana (35) 2 and watermelon (35) 20 . It seems that some aquaporin genes physically cluster very close to each other in certain regions, suggesting the occurrence of more than one gene tandem duplication event in the evolutionary history of melon, as was previously studied in tobacco plants 58 . There is an accumulation of the PIP2 genes in chromosome 5. In the same way, the co-location of CmNIP2;1 and CmNIP2;2 in chromosome 4 and of CmNIP5;1 and CmNIP5;2 in chromosome 9 could allude also to genomic duplication events. Despite this, further investigations are needed to clarify the evolutionary events that have occurred across the aquaporin family in C. melo.
The protein features, such as the mean Mw, of each subfamily also showed similarities to those of A. thaliana 59 . However, a notable difference is shown in the SIPs subfamily, the Mw being greater in melon than in A. thaliana; this could be due to the larger sequence found in CmSIP2;1, compared with AtSIP2;1. Furthermore, when the protein sequence of CmSIP2;1 was compared with those of homologous SIP2;1 proteins in other species (watermelon, arabidopsis and cucumber) 2,20,47 , all the sequences aligned from the second methionine residue of CmSIP2;1. This suggests that this site is the starting point of translation, resulting in a protein of 238 aa and 25.99 kDa, thus equating the results to those obtained in arabidopsis for this subgroup. However, this hypothesis needs to be studied in order to clarify all the protein features previously predicted. In addition, TIPs were not found to be smaller than PIPs and NIPs, but most of them were more acidic (with lower Pi) than PIPs, NIPs, SIPs and XIP. The Pi could reflect a functional constraint imposed on MIPs 60 . Indeed, sequence analysis in arabidopsis revealed that the cause of this difference in Pi lies in the C-terminal regions, which are more basic in PIPs and NIPs than in TIPs. Therefore, it is possible that phosphorylation sites or sorting signals in the C-terminal regions form part of the hypothesised functional constraint on the sequences 60 . www.nature.com/scientificreports/ The aquaporins monomers found a highly conserved structure, with six transmembrane helices in all CmAQPs (Table 1, Fig. 3). These strongly conserved regions has been reported to be likely constrained to maintain the structural integrity of the aquaporin monomer and the conservation of critical residues which was essential for tetramer formation 61 . In this way the subcellular location of aquaporins, PIPs, NIPs and XIPs are usually located in the plasma membrane, while TIPs and SIPs are normally localised to the tonoplast and endoplasmic reticulum (ER), respectively 62,63 . However, our results showed that the PIPs, NIPs, SIPs and XIP were located principally in the plasma membrane (Table 2). This was expected, with the exception of the SIPs-that were not located in the ER by any of the programs used, although this is their most probable location. Nevertheless, a few of NIP proteins were located not only in the plasma membrane but also in the tonoplast and the TIPs showed a diverse range of possible subcellular locations depending on the chosen prediction software (chloroplast, mitochondria, plasma membrane and cytoplasm). Significantly, CmTIP5;1 was the only aquaporin assigned to chloroplasts, while CmTIP3;1 was localized in mitochondria ( Table 2). As cell compartmentalization has represented the main driving force in the diversification of aquaporins in plants, no location should be excluded without specific analysis and more studies must be made to confirm these preliminary data. Nevertheless, NIPs isoforms have not yet been found in the tonoplast, so this does not appear to be their most likely location. A role for AtTIP5;1 as a nitrogen transporter in mitochondria has been proposed 64 and different isoforms of PIPs and TIPs have been detected in chloroplasts 65,66 ; a role in the transport of CO 2 , water and H 2 O 2 in chloroplast membranes could be beneficial to the plant dynamics 62 . So, the predicted localization of CmTIP3;1 in mitochondria and of CmTIP5;1 in chloroplasts cannot be rejected, even if this location for CmTIP5;1 does not seem to coincide with the gene expression data (Fig. 4c), which place it mostly in roots. (Table 3). The function of PIPs in water transport is very specialised and their structure is highly preserved 67 . All PIPs have the same NPA motifs and ar/R selectivity filters residues (F, H, T and R in H2, H5, LE1 and LE2, respectively) ( Table 3), directly related to their main role in water transport 7 , while phylogenetic and functional studies suggest that they are equally capable of transporting CO 2 10,12 . The FPs are also very similar in most PIP aquaporins from both subgroups (PIP1s and PIP2s). The S-A-F-W residues are well conserved in all cases in positions P2 to P5 respectively. In nine of the 12 PIPs, P1 was a Q (Gln) residue; the exceptions were CmPIP1;2 with E (Glu), CmPIP2;9 with K (Lys) and CmPIP2;10 with M (Met) ( Table 3). The comparison of melon FPs with homologues in other species predicted the possible transport of H 2 O 2 by almost all melon PIPs 12 . The only exception was the unusual CmPIP2;9, whose aa residues in the ar/R filter (F-N-A-R) and in the P1-FPs (K) differed from those of the PIP2s characterised to date and whose particular sequence has not been previously described. The sequence homology analysis (Supplementary Table S1) of PIPs that have been shown to transport B, urea and As resulted in the prediction of matching transport by orthologues in C. melo 11 . According to this, CmPIP1s could be able to transport boric acid and/or urea, like their orthologues ZmPIP1;1 and ZmPIP1;5 31,32 , while CmPIP2;10 could transport As, as OsPIP2;6 does 45 (Supplementary Figure S2).

Analysis of possible functions in solutes transport
Based on the ar/R filter 68 , the CmTIPs were classified in four groups attending to their homology. Group I is formed by CmTIP1;1, CmTIP1;2 and CmTIP1;3, Group IIa is formed by CmTIP2;1 and CmTIP2;2, Group IIb is constituted by CmTIP3;1 and CmTIP4;1 and Group III has only one member, CmTIP5;1. In general, TIPs seem to have developed the capacity to transport nitrogenous compounds 69 . The great importance of the H2 and H5 positions (H-I) with a non-polar LE1 (A/G) in the groups I, IIa and IIb has been related with NH 3 transport 30 . Many TIPs have been characterised as urea transporters 32,33,46,70 . In the TIPs group I, the substitution of the typical R by V in the LE2 position has been proved to be involved in the transport of not only NH 3 and but also H 2 O 2 , in mutagenic studies 30,71,72 . Regarding the FPs, the CmTIP1s presented slight differences, and a direct homology with residues of the other plants studied was not found (Supplementary Table S1). However, the phylogenetic analysis clearly related the CmTIP1s with the TIP1s of A. thaliana (Fig. 1). All the AtTIP1s were predicted to transport urea and H 2 O 2 10-12 and this has been proven experimentally in heterologous systems 33,36,70 , while NH 3 transport has only been tested in mutagenic studies 72 . All this suggests that C. melo TIP1s are able to transport these solutes. The clade including CmTIP1;3, AtTIP1;3, ZmTIP1;2 and OsTIP1;2 diverged earlier, followed by the separation of CmTIP1;2 that presents a greater evolutionary divergence from AtTIP1s but is closer to OsTIP1;1 and ZmTIP1:1 (Supplementary Figure S2). Interestingly, it has been shown that OsTIP1;2 can transport glycerol 73 and that ZmTIP1;1 and ZmTIP1;2 can transport boric acid, in addition to NH 3 , urea and H 2 O 2 46 . So, CmTIPs1 might be able to transport B and CmTIP1;3 might be able to transport Gly also (Table 3). Regarding Group IIa (TIP2s), their ar/R region has been directly related to the transport of NH 3 and urea 28,33,72,74 . A specific amino acid, histidine (H), in loop C seems to be key for the de-protonation of NH 4 + , which would allow the transport of NH 3 independent of pH 74 . This residue is present in CmTIP2s and CmTIP4.1, but not in the other TIPs subgroups, which lack the ability to extract a proton from an NH 4 + ion; thus, the passage of ammonia relies on the pH-dependent concentration of uncharged NH 3 in the medium 71 . The transport predictions for FPs indicate that they can transport not only NH 3 and urea but also H 2 O 2 10-12 ; their relationship with AtTIP2s ( Fig. 1) also suggests this capability, while the ability of orthologues to transport H 2 O 2 has only been tested in assays with mutants 75 . Phylogenetic analysis of the TIPs group IIb showed that CmTIP3;1 is closely related to AtTIP3s (Fig. 1), although there is a substitution in FPs P3 (A is replaced by S) ( Table 3). The ar/R analysis supports the NH 3 and urea transport capacity 28,30,72 . As stated previously, CmTIP4;1 and CmTIP2s share all FPs, the His in loop C 74 and the same ar/R region (with the only exception of G instead of A in LE1) (Table 3), pointing to selective NH 3 transport. Urea transport by CmTIP4;1 is supported by mutant studies 28 and has been shown to be transported by the orthologous AtTIP4.1 33 . The prediction based on the phylogenetic framework, and including FPs, also indicated the ability to transport H 2 O 2 12 . CmTIP5;1 is the only aquaporin in group III and its orthologue AtTIP5;1 has been shown to transport urea 70 . Their FPs differ only in the P1 position (Table 3), which in C. melo Scientific Reports | (2020) 10:22240 | https://doi.org/10.1038/s41598-020-79250-w www.nature.com/scientificreports/ is a residue of I instead of V (both non-polar); therefore, the transport of urea by CmTIP5;1 seems to be possible and is supported by mutagenic analysis of the ar/R region 28 . The NIPs form a monophyletic group 76 and overwhelming evidence supports a main role for NIPs in metalloids transport 77 beside glycerol-including B, Si, As and Sb [reviewed in Pommerrenig et al. (2015)]. Interestingly, all these compounds, in their uncharged states, have a geometry that is substantially similar to a conformation of glycerol in a retracted state; this would readily allow adaptation of the pore of the original aquaglyceroporin NIPs to the transport of these compounds 78 . Eight NIPs have been found in C. melo and have been classified into the three subgroups of NIPs based on their ar/R filter 79 . CmNIP1;1 and CmNIP4;1 belong to Group I. CmNIP5;1, CmNIP5;2, CmNIP6;1 and CmNIP7;1 belong to group II, although some of their compositions vary from the typical ar/R residues (discussed below), and CmNIP2;1 and CmNIP2;2 belong to Group III. The typical W-V-A-R in H2, H5, LE1 and LE2, respectively, are the characteristic ar/R residues in Group I and have been directly related to water and glycerol transport 29 . Both CmNIP1;1 and CmNIP4;1 have these characteristic motifs ( Table 3). The phylogenetic tree places both in the same clade along with AtNIP1;1, AtNIP1;2, AtNIP2;1 and AtNIP3;1 (all in a branch with CmNIP1;1), AtNIP4;1 and AtNIP4;2 (both in the same branch as CmNIP4.1) (Fig. 1). Experimental studies showed that AtNIP1s were able to transport As, Sb, H 2 O 2 and Gly 37,42,75,80 . CmNIP1;1 appears to have separated from AtNIP1s earlier, in an evolutionary branch that it shares with maize and rice NIP1s (Supplementary Figure S2). OsNIP1;1 can transport As 38 while ZmNIP1;1 has been shown, in heterologous systems, to transport H 2 O 2 and glycerol; unexpectedly, it is also able to transport boric acid efficiently 46 . Thus, CmNIP1;1 might transport As, Sb, H 2 O 2 and Gly, while B cannot be ruled out. Interestingly, while CmNIP4;1 belongs to group I and has the same ar/R motifs, its FPs are more similar to those of NIPs group II. This increases the possibility that it transports other interesting solutes, such as those that can be transported by NIP2s (discussed below). Despite this, the homology-based prediction did not allow a specific transport role to be assigned to CmNIP4;1, beyond glycerol transport compatibility 29 , due to its clear differences in FPs from its orthologues (Supplementary Figure S2 & Supplementary Table S1). Strong evidence supports the role of the NIPs II group in boric acid transport 34,41,81 . Among the C. melo NIPs belonging to this group, the most important variation is in H2; in CmNIP5;2 it is an S residue and in CmNIP6;1 a T residue (Table 3). However, both CmNIP5;2 and CmNIP6;1 conserve FPs identical to those of CmNIP5;1 and homologous to those of other aquaporins of group II found in arabidopsis, rice and maize (Supplementary Table S1), all of them belonging to the same clade in the phylogenetic tree (Supplementary Figure S2). CmNIP5;1 and CmNIP5;2 are very close to AtNIP5;1, which has been proved to transport As, Sb and B in experimental assays 26,34,82 . CmNIP6;1 is related directly to AtNIP6;1, the latter being able to transport not only As, Sb and B but also urea and Gly 27,41,82 (Fig. 1). The prediction based on the phylogenetic framework, and including FPs, assigned to CmNIP5;1 the ability to transport As, Sb and B, as AtNIP5;1 does. This As and B transport is supported by mutagenic analysis 26 and its ability to transport Gly and/or urea must also be considered 27 . The variation in H2 of CmNIP5;2 (Table 3) prevented the prediction of any transport ability, although it should not be ruled out that it fulfils functions similar to those of CmNIP5;1. The same applies to CmNIP6;1, for which the prediction points only to Sb transport, although this is most likely not its only function. Regarding CmNIP7;1, its phylogenetic divergence does not separate it from group II of the NIPs and its proximity to AtNIP7;1 suggests that it can perform similar functions (Fig. 1). AtNIP7;1 can transport As, Sb, B, urea and Gly 73,82 . The possible transport of Gly and/or urea is also supported by studies with point mutations 27 and by the presence of a specific Y64 residue (Y81 in AtNIP7;1) which has been directly related to gating and regulation of urea and Gly transport 83 . The aquaporins CmNIP2;1 and CmNIP2;2 are included in the NIP III group 79 . They possess a unique ar/R selectivity filter due to the small size of the aa in the H2 and H5 positions, giving the largest pore diameter described in aquaporins, that may allow the passage of very large solutes such as Si (4.38 Å) 84 . Indeed, the transport of Si has been specifically associated with the motifs that are present in CmNIP2;1 85 . Si transporters are absent from the Brassicaceae and, therefore, NIP III aquaporins have not been found in A. thaliana, but they are present in all silicon-accumulator plants. Phylogenetic analysis assigned CmNIP2;1 and CmNIP2;2 to the same clade as the Z. mays and O. sativa NIP2s (Supplementary Figure S2), strongly suggesting that all members of the C. melo NIPs III group are Si transporters. In addition to silicon, members of the NIPs group III from Z. mays and O. sativa have been proved experimentally to transport As, Sb, B and urea 39,46,79,[86][87][88] ; also, ZmNIP2:1 was able to transport small amounts of H 2 O 2 and Gly in heterologous expression experiments 46 . Thus, according to the analysis of orthologues, the transport of most of these solutes could occur via C. melo NIP2s and the prediction based on homologies also points to such transport by CmNIP2;1 [10][11][12] . The case of CmNIP2;2 is peculiar since it clearly belongs to the same protein family but has a variation in the H2 position of the ar/R selectivity filter (C instead of G), and also its FPs are slightly different (S in P2, instead of T) from those of the NIPs III group (Table 3); hence, a prediction could not be performed based on homologous genes. The substitution of G in H2 by another non-polar small residue (A) did not affect the ability of OsNIP2;1 to transport Si, As or B 26 , so it is reasonable to propose that CmNIP2;2 can also transport some of these solutes.
Very little is known about the function of SIPs in planta. They are localised in the ER and facilitate water transport 89 but no data on transport of non-aqueous substrates are available. Despite this, homologues-and structure-based analyses predict that the SIP1s of A. thaliana, Z. mays and O. sativa are capable of transporting urea, while the SIP2s could transport H 2 O 2 and As 12 . Phylogenetic analysis showed that the SIPs of C. melo are clearly related to their orthologues in A. thaliana (Fig. 1), although the variation that both CmSIP1;1 and CmSIP2;1 show in H2 impairs the direct assignation of solutes transport to these proteins ( Table 3).
The absence of XIPs in the Brassicaceae and monocots makes it impossible to compare the sequence found in C. melo with those of A. thaliana, Z. mays and O. sativa. For this reason, we compared the important residues with those of XIPs in other species 90 . CmXIP1;1 belongs clearly to the XIP-A subfamily, which does not yet have known transport functions. The main particularity of Cucumis (both sativus and melo) XIPs is that they present an SPI/SPA motif rather than a NPA/NPA motif, while conserving other MIP-specific sequence features. It is Scientific Reports | (2020) 10:22240 | https://doi.org/10.1038/s41598-020-79250-w www.nature.com/scientificreports/ usual to find variations in the first NPA motif of XIP1s-as they have been observed in other species, such as Ricinus communis L., Lotus japonicus L., Prunus persica L. and Glycine max L. 90 -but it is quite odd to find a mutation in the second NPA motif. Also, in the first NPA motif, the A is substituted by I, which is also quite rare and has been observed only in Cucumis to date. It is still unknown what kind of functional effect this residue change might have on the protein.

Constitutive expression of CmAQPs and possible implication in melon physiology. Our results
clearly show that the main aquaporin expressed in C. melo was CmTIP1;1, by a wide margin with respect to the other isoforms (Fig. 4a). The water transport capacity has been maintained in TIPs, especially in the TIP1s, for which the magnitude stands out, being mostly responsible for the high permeability of the tonoplast 91 . CmTIP1;1, therefore, should be mainly responsible for the water transport between the cytoplasm and vacuole in melon plants. The other TIP1s in C. melo showed very slight expression (Fig. 4c). If this is due to a redundancy of functions because of the importance of CmTIP1;1, or to their specialisation in other tissues or environmental conditions, remains to be investigated. Just below CmTIP1;1, the PIP1s isoforms were the aquaporin genes most expressed in our studies, in both roots and leaves (Fig. 4a). Other studies showed patterns similar to ours. For example, in maize, in the central zone of the leaf, the aquaporin gene most expressed was shown to be ZmPIP1;1, with PIP1s contributing more than 70% of the expression of PIPs in leaves 92 . The characteristic structure of PIPs is also clearly compatible with CO 2 and H 2 O 2 transport besides that of water. So, certain PIP1 isoforms might specialise in such functions, depending on their location, and some PIP1s were shown to be directly related with CO 2 transport in mesophyll cells and chloroplasts 66,92,93 . This allows us to propose a function of CmPIP1s in the transport of both water and CO 2 in melon, which could explain their extremely high expression in leaves. The PIP2s isoforms, discounting their specific function as water transporters, have been related with ROS signalling through their regulation by and transport of H 2 O 2 36,75,94 . Attending to the expression profile, CmPIP2;6 seems to have also a fundamental role in leaves (Fig. 4b), as occurs in other species such as Arabidopsis, in which one PIP2 (AtPIP2;6) is predominant in leaves 95 , pointing to a possible implication in H 2 O 2 or even CO 2 transport besides that of water. Indeed, it has been shown that, despite the high homology among the different PIPs, some isoforms of plants possess the ability to transport water but not CO 2 (NtPIP2;1) or vice versa (PIP1, NtAQP1), the proportion of each of them in the tetramer allowing the passage of water or CO 2 to be regulated. This opens the door to the involvement of the fifth, central pore in the transport of gases and to possible regulation of the transport by competition between subunits in the formation of the tetramer 5 . How tetramers can be established preferentially with one concrete isoform, to regulate such transport, could be of interest in future studies to see if CmPIP2;6 is preferentially selected in tetramers with PIP1s involved in CO 2 conduction.
In addition, PIP1s also had high expression in roots along with CmPIP2;10, which is probably the main water transporter in C. melo, especially in roots, followed by CmPIP2;2, CmPIP2;3 and CmPIP2;4 (Fig. 4b,c). As we show in this study, and in general under optimal conditions, there is greater expression of aquaporins in the roots 51,96,97 , since these are the point of entry of water and nutrients into the plant, and the predominance of a specific PIP2 isoform also seems to be habitual. As an example, ZmPIP2;5 is the main PIP gene expressed in the roots of maize and has been shown to be essential for water to cross the Casparian barrier in the exodermis 51 . The heterotetrameric conformation that includes PIP2 and PIP1 subunits has been shown to improve water transport, versus heterotetramers consisting exclusively of PIP2s 98,99 . This could be a valuable source of water transport regulation in the roots of melon plants as the interaction of the two isoforms is essential for PIP1s to access their active location in the plasma membrane 99 . Interestingly, the orthologues of CmPIP1;1, CmPIP1;2 and CmPIP2;10, the most expressed PIP isoforms in melon roots, have been found to transport other solutes such as urea, B and As 31,32,45 . The occurrence of this transport through these isoforms in melon could be of great interest to understand their expression patterns. Apart from CmTIP1;1, within the TIPs, only CmTIP3;1 seems to have a significant presence under optimal conditions, in both roots and leaves (Fig. 4b). Given that the prediction of the subcellular location (Table 2) places CmTIP3;1 in the mitochondria and the well-known ability of TIPs to transport nitrogen compounds then, if this location is confirmed, they could also play a role in the transport of NH 3 produced by photorespiration 63 . Attending to the CmNIPs, the expression of CmNIP5;1 and CmNIP5;2 was similar, being low but significant in both roots and leaves (Fig. 4b). Strong evidence supports the role of the NIPs II group in boric acid transport 34,41,81 and this is probably the main function of the CmNIP5s. The members of the NIPs III group have been defined clearly as Si transporters 100 and the CmNIP2s were the most expressed NIPs genes in melon plants, especially CmNIP2;1, whose expression was higher in leaf tissues, while CmNIP2;2 was expressed mostly in roots (Fig. 4b). In rice it has been shown that OsNIP2;1 is the main Si transporter involved in Si uptake, while OsNIP2;2 is mostly implicated in the unloading of Si from xylem vessels in leaves 35,43 . Similarly, the stronger presence of CmNIP2;1 in leaves suggests a Si-unloading function in leaves, while the two isoforms had the same expression in roots, implicating both in the uptake and translocation of Si in roots. However, although neither B transport by CmNIP5;2 nor Si transport by CmNIP2;2 could be directly predicted based on previous analyses (Table 3), these are probably their main functions. This idea is supported by the tight tandem that each of them forms with its respective paralog in chromosomes 9 and 5, respectively. These probably represent the result of gene duplication events that yielded differences in discrete functions or sublocalisations within the same tissue, or even between different tissues, as seems to have been the case in the NIP2s.
Finally, it is worth noting the expression of the aquaporin CmSIP1;1 in roots, which was at the level of many PIPs and more than half of the isoforms analysed. Thus, although specific functions are not known, it seems to have some importance in the physiology of melon roots.
Comparison of expression profiles in RNA-seq versus qPCR. Quantitative reverse transcription PCR (RT-qPCR) was the tool used most widely for the quantification of gene expression due to its sensitivity and Scientific Reports | (2020) 10:22240 | https://doi.org/10.1038/s41598-020-79250-w www.nature.com/scientificreports/ precision 101 , being considered the most important medium-performance gene-expression analysis technology 102 , and it should be used to validate RNA sequencing, which is the most widely used tool today. Based on this, we decided to compare the expression of the 31 aquaporins genes of C. melo var. Grand Riado with the RNA-seq data available in the databases, which refer to another variety, in this case Cantaloupe (melon) 56 . We found numerous similarities between our qPCR analyses and the RNA-seq data. Among them, the aquaporins with higher levels of expression showed better correlations for one or both tissues; this was the case for PIP1;2, PIP2;1, PIP2;6, PIP2;10, TIP1;1, TIP4;1, NIP2;1 and NIP5;1. The differences in expression between the qPCRs and RNA-seq can be explained, since the sensitivity of the qPCRs decreases from cycle 34 onwards and most of the differences were found in aquaporins that appeared from this cycle onwards due to their low expression. Other mismatches between the two analyses can be explained by the use of different melon varieties. This could have led to changes in the expression pattern of some aquaporins-as seen in other plant varieties, where differences have been found, even between different cultivars 103 . Lastly, the growth conditions were different: Latrasse et al. (2017) used constant humidity (60%) and a temperature of 21-27 ºC in a glasshouse (for melon plants grown on rock compost) and also grew plants in the field (in soil), while we used hydroponic culture, at 60-80% humidity and a temperature of 20-25 ºC 56 .
Despite the discrepancies, a clear relationship is seen regarding the importance of certain isoforms in C. melo. Clearly, CmTIP1;1 is revealed as the most important osmoregulator in the tonoplast under optimal conditions. The PIP1s are of greater importance, in both roots and leaves, CmPIP1;1 being the predominant isoform. CmPIP2;6 is very important in leaf tissues, while in the case of the variety Cantaloupe its function seems to be shared with CmPIP2;10, which had a strong presence in leaves. In roots, the Cantaloupe aquaporins isoforms that seem to have the major role in water uptake are CmPIP2;4, followed by PIP2;10, PIP2;2 and, finally, PIP2;1. Of the NIPs, NIP5;1 and the NIP2s are undoubtedly important in melon plants. In the variety Cantaloupe, the presence of NIP5;1 is extremely striking in roots, being one of the most expressed isoforms, while in Grand Riado both NIP5 isoforms may share redundant functions. The genetic duplication of the NIP2s in a narrow chromosomal tandem and their matching transport abilities, together with their differentiation with respect to the expression in roots and leaves, suggest that they have acquired tissue-specific roles, these possibly being different in each of the two varieties. Finally, both analyses indicated that SIP1;1 plays a constitutive and important, but still unknown, role. It seems that, in each variety, slight changes in the weight of different paralogs have developed and these changes could be related to the transport function and/or localisation. The abundant solute transport possibilities for each aquaporin isoform could have played an important role in the specialisation within each variety, depending on the needs of each of these two varieties throughout their evolution and development in different environments.

Concluding remarks
In this work, the sequences of aquaporins genes have been assessed in melon. The comparative analysis with related plant species was very useful to predict the possible roles of some of the aquaporins in the transport of certain solutes. This analysis together with qPCR carried out in roots and leaves revealed the constitutive expression of almost all aquaporins genes in both organs. However, the role of each aquaporin must be elucidated in future work. Therefore, the exhaustive design of all the primers presented here is intended to serve as the basis for future research into all aquaporins in melon, regarding their roles in plant development and in the response to stress conditions. Finally, this work will serve as a basis for other researchers who wish to carry out a comprehensive analysis of the genome, explaining and clarifying a simple way to carry out this type of research.

Data availability
The seeds were kindly provided by SAKATA SEED IBERICA, S.L.U. and all data generated or analysed during this study are included in this published article and its supplementary information files. Original data from qPCR are available from the corresponding author on reasonable request.