Freshly harvested vegetative and reproductive organs/tissues (1 month-old plantlets; pool of leaves of different ages; roots from 1 and 6 month-old plantlets; stems from 1 and 6 month-old plantlets; vascular cambium samples from 6 year-old trees; flower; flower buds; and fruits) were obtained from Eucalyptus grandis essentially as described previously (Rodrigues et al., 2013). After harvesting, fresh samples were frozen immediately in liquid nitrogen until total RNA extraction. Two-month-old plantlets of a commercial clone of E. grandis, kindly provided by Suzano Papel e Celulose SA, Brazil, were used in the osmotic stress assays.

To identify members of the TIP gene subfamily in Eucalyptus, annotated TIP sequences from Arabidopsis (Johanson et al., 2001; Alexandersson et al., 2005) were used as queries in BLAST searches against the E. grandis BRASUZ1 genome assembly (v2.0) available at Phytozome¹. To refine our analysis, additional BLAST searches were performed using previously identified poplar TIP sequences (Gupta and Sankararamakrishnan, 2009) as drivers. The identified gene products were subsequently validated using the Pfam domain annotation (Finn et al., 2016). The exon/intron organization of the mined EgTIP genes was determined using the Eucalyptus gene models annotated in Phytozome and the Gene Structure Display Server² (GSDS v2.0). The existence of duplication of the EgTIP genes was investigated using the locus-search module of the Plant Genome Duplication Database³ (PGDD). The ratio of non-synonymous (Ka) to synonymous (Ks) substitution rates of evolution were also calculated using PGDD. Searches for the presence of putative cis-regulatory elements within the EgTIP promoter regions were carried out using the Plant Care and PLACE databases. Predictions of the subcellular location of the EgTIP proteins were performed using Plant-mPLoc⁴ (Chou and Shen, 2010).

The deduced aa sequences of the identified EgTIPs were aligned with known TIPs from A. thaliana, P. trichocarpa, Quercus petraea, and V. vinifera using CLUSTALX2⁵. Phylogenetic relationships were inferred using the neighbor-joining method with 1000 bootstrap replicates implemented by the MEGA 7 software package⁶. Branches with <50% bootstrap support were collapsed. Accession numbers of the sequences used in the phylogenetic analyses are as follows: AtTIP1.1 (At2g36830), AtTIP1.2 (At3g26520), AtTIP1.3 (At4g01470), AtTIP2.1 (At3g16240), AtTIP2.2 (At4g17340), AtTIP2.3 (At5g47450), AtTIP3.1 (At1g73190), AtTIP3.2 (At1g17810), AtTIP4.1 (At2g25810), AtTIP5.1 (At3g47440), QqTIP1.1 (JQ846274), QqTIP2.1 (JQ846275), QqTIP2.2 (JQ846276), VvTIP1.1 (GSVIVP 00018548001), VvTIP1.2 (GSVIVP00000605001), VvTIP1.3 (GSVIVP00022146001), VvTIP1.4 (GSVIVP00024394001), VvTIP2.1 (GSVIVP00034350001), VvTIP2.2 (GSVIVP00012703001), VvTIP3.1 (GSVIVP00013854001), VvTIP4.1 (GSVIVP00032441001), VvTIP5.1 (GSVIVP00029946001), VvTIP5.2 (GSVIVP00019170001), PtTIP1.1 (549212), PtTIP1.2 (833283), PtTIP1.3 (822504), PtTIP1.4 (656044), PtTIP1.5 (667870), PtTIP1.6 (589502), PtTIP1.7 (558321), PtTIP1.8 (828458), PtTIP2.1 (548890), PtTIP2.2 (645978), PtTIP2.3 (817166), PtTIP2.4 (676397), PtTIP3.1 (584517), PtTIP3.2 (811826), PtTIP4.1 (561759), PtTIP5.1 (414059), PtTIP5.2 (423803). The Phytozome IDs of the mined E. grandis TIPs are: EgTIP1.1 (Eucgr.K02438), EgTIP1.2 (Eucgr.J00051), EgTIP1.3 (Eucgr.B02403), EgTIP1.4 (Eucgr.J02074) EgTIP2.1 (Eucgr.F03054), EgTIP2.2 (Eucgr.D02090), EgTIP2.3 (Eucgr.D02507), EgTIP3.1 (Eucgr.K02339), EgTIP3.2 (Eucgr.H00668), EgTIP4.1 (Eucgr.C02914), and EgTIP5.1 (Eucgr.K00164).

The deduced aa sequences of EgTIPs representing each TIP group (EgTIP1.1; EgTIP2.1; EgTIP3.1; EgTIP4.1, and EgTIP5.1) were employed for the construction of the models. These structures were predicted based on the alignment data generated by the program Phyre2 that uses homology detection methods to build 3D models (Kelley et al., 2015). Considering the quality of the structure, the crystallographic structure of the A. thaliana aquaporin TIP2.1 (AtTIP2.1) at 1.18 Å resolution (PDB ID: 5I32, Chain A) (Kirscht et al., 2016) was chosen as template. The program MODELLER 9v16 (Sali and Blundell, 1993; Marti-Renom et al., 2000) was used to generate the protein models and rank them according to the DOPE (Discrete Optimized Protein Energy) score. The best TIP models were selected based on the stereochemical parameters using the program RAMPAGE (Lovell et al., 2003). The ConSurf Server (Landau et al., 2005; Ashkenazy et al., 2010) was used for conservation analysis and all images were generated using the PyMOL program (Schrödinger, 2011).

The Eucalyptus plantlets (two months old) were transferred to hydroponic culture by means of a floating system and maintained in an aerated Hoagland’s nutrient solution (75%) at pH 6.0 (osmotic potential of -0.1 MPa) before stress imposition. After an acclimatization period, 50% of the tested plantlets were stressed by adding 215 g of polyethylene glycol 8000 (PEG; Sigma, USA) to 1 L of culture medium to induce osmotic stress, while the remaining were maintained in Hoagland’s nutrient solution (75%) throughout the assay as control. The osmotic potential of the PEG 8000-treated solution was -0.6 MPa (WP4-T, Decagon Devices, Inc., England). A compressed air system was used to homogenize the solutions in both the tray and the PEG container to avoid anoxia of the roots.

Both control and stress treatments comprised 15 plantlets per tray, and three trays per treatment (control/treated plantlets), representing three biological replications. To carry out the gene expression analyses, the roots, stems and leaves from three randomized plantlets per tray were collected after 6 h (short-term), 12 h (medium-term), and 24 h (long-term) of PEG-treatment. Organ samples from untreated control plantlets were simultaneously collected at each time-point. To reduce plant-to-plant variation, each collected group of organ samples was pooled before RNA extraction.

The relative expression of the EgTIP genes was assessed using quantitative real-time RT-PCR (RT-qPCR). Total RNA extraction and cDNA synthesis were performed as previously described (Rodrigues et al., 2013). The qPCR analyses were carried out using Power SYBR Green Master Mix (Applied Biosystems) and a StepOnePlus Real Time PCR System (Applied Biosystems). The cycling conditions were as follows: 5 min at 95°C, followed by 45 cycles of 15 s at 95°C and 60 s at 60°C. Each reaction was performed in triplicate in a total volume of 10 μl, and contained 20 ng of cDNA and 0.2 μM of each EgTIP-specific primer (Supplementary Table S1). Among five selected E. grandis reference genes (Actin, GAPDH, Cdk8, Transcription elongation factor s-II, and Aspartyl-tRNA synthetase; de Oliveira et al., 2012) tested, actin was the most stable and thus employed as endogenous control (Supplementary Table S1; Goicoechea et al., 2005). Cycle threshold (Ct) values were obtained for each sample, and relative quantification was determined using the 2^-ΔΔCt method as described (Livak and Schmittgen, 2001) Amplification efficiencies were derived from the amplification plots using the program LinRegPCR (Ramakers et al., 2003). A value of two was used in calculations. For the expression analysis under osmotic stress, the relative expression of each EgTIP at a given time point was determined as the fold change of its expression under treated condition relative to its expression under control condition.

Relative expression data were analyzed using the Relative expression software tool (REST 2009) and differences with p-values <0.05 were considered statistically significant. The heatmaps were created using the function heatmap.2 from plot package at R environment. For the organ/tissue-specific gene expression analyses, normalization was performed using 1-month-old plantlets as control sample as previously reported (Yu et al., 2014). According to Yu et al. (2014), these plantlets represent a highly stable and less variable sample that contain the main investigated organs/tissues. Additional expression analyses were performed using publicly available RNA-Seq data⁷ generated for six different vegetative organs/tissues (young and mature leaves, shoot tips, phloem, immature xylem and xylem) from a 6-year-old E. grandis X E. urophylla hybrid clone (Mizrachi et al., 2010). In this case, transcript abundances were expressed as units of normalized FPKM (fragments per kilobase of exon per million fragments mapped), that were used for heatmap construction.

BLAST searches in the E. grandis genome identified 11 putative TIP coding sequences (Table 1), thus indicating that the Eucalyptus TIP subfamily has 11 members (refereed as EgTIP) and not 5 as previously suggested (Rodrigues et al., 2013). To analyze in more detail, the evolutionary relationships between the identified EgTIPs and those present in other plant species, an unrooted phylogenetic tree was constructed using the predicted protein sequences (Figure 1). Based on their phylogenetic relationships, the EgTIPs could be individually assigned to each of the classical TIP1–5 groups (Johanson et al., 2001), and were named accordingly (EgTIP1.1, EgTIP1.2, EgTIP1.3 and EgTIP1.4; EgTIP2.1, EgTIP2.2 and EgTIP2.3; EgTIP3.1 and EgTIP3.2; EgTIP4.1; and EgTIP5.1). In general, the corresponding groups, with the exception of group 1, were supported by moderate to high bootstrap values. The identity at the aa level varied according to the TIP group considered: 79–88% identity between EgTIP1s, 76%–47% identity between EgTIP2s and 77% between EgTIP3s. Taking into account the tree topology, EgTIP2.3 was grouped with known TIP5 isoforms, being closely related to VvTIP5.2 (61% identity at the aa level). Nevertheless, the inspection of the ar/R selectivity residues of both VvTIP5.2 and EgTIP2.3 (see below) revealed a closer relationship with TIP2 isoforms. Therefore, based on its selectivity filter, EgTIP2.3 was assigned here as a TIP2 as originally proposed in Phytozome instead of a TIP5. In fact, the small subgroup formed by VvTIP5.2 and EgTIP2.3 is branched between the groups encompassing TIP5.1 and TIP2 isoforms. In this regard, EgTIP2.3 shares 44 and 47% aa sequence identity with EgTIP5.1 and EgTIP2.2, respectively. As observed in other plant species, the majority of the identified EgTIPs belongs to group 1 (four members) followed by group 2, 3, and 5 with two members each. The phylogenetic tree also revealed the presence of closely related gene pairs (EgTIP1.3/EgTIP1.4 and EgTIP3.1/EgTIP3.2, for example), suggesting possible gene duplication events.

Analysis of the genomic structure of the annotated EgTIP genes revealed a two-intron/three-exon organization, except for EgTIP1.1 that contained two exons disrupted by a single intron (Figure 2). These introns showed variations in length and position. Consistent with this, most members of the TIP subfamily of dicot species share such a well-conserved two-intron/three-exon pattern (Regon et al., 2014). Regarding their chromosomal distribution, the EgTIP genes are located on seven (2, 3, 4, 6, 8, 10, and 11) out of the eleven eucalyptus chromosomes (Table 1). In general, one or two genes are found per chromosome. In this context, EgTIP1.1 and EgTIP4.1 (on 3), EgTIP1.2 and EgTIP1.4 (on 10), EgTIP2.2 and EgTIP2.3 (on 4) and EgTIP3.1 and EgTIP5.1 (on 11) are among those located on the same chromosome.

An additional locus search at the PGDD website was performed to evaluate possible EgTIP gene duplication. According to the resulting data, the following gene pairs were observed: EgTIP1.3/EgTIP1.4, EgTIP2.1/EgTIP4, EgTIP3.1/EgTIP3.2, and EgTIP2.2/EgTIP2.3. Among them, only EgTIP2.2 and EgTIP2.3 are located on the same chromosome (4), and may possibly represent tandemly duplicated genes. To investigate the mechanisms involved in gene duplication evolution after divergence, we calculated the Ka/Ks ratio of each EgTIP gene pairs. The results (Supplementary Table S2) revealed that the Ka/Ks ratios of all four EgTIP gene pairs were less than 0.4, implying that these genes had undergone purifying selection pressure with limited functional divergence after duplication.

The deduced EgTIP proteins ranged in size from 243 (EgTIP4.1) to 262 (EgTIP3.1) aa residues, with molecular weights (MW) and theoretical isoelectric point (pI) values raging between 25 and 27.8 kDa and 4.68 and 6.81, respectively (Table 1). These proteins are all predicted to contain six transmembrane α-helices (H1–H6) and five loops. Concerning their sub-cellular localization, all EgTIPs are predicted to be located in the vacuole, with the exception of EgTIP5.1 and EgTIP2.3, that are also localized at the cell membrane (Table 1).

To evaluate structural similarities and differences, theoretical models for EgTIPs representing each classical TIP group were constructed using a comparative modeling approach and the crystallography structure of AtTIP2.1 at 1.18Å resolution (Kirscht et al., 2016) as a template. The sequences shared 53–85% aa identity and 90–95% of coverage. As expected, all models contained the six α-helices and interhelical loop regions as seen throughout the aquaporin family. The best models showed good stereochemical quality with 96.9% of the residues in favored regions, 2.6–1.7% and 0–0.9% (corresponding to two residues) in allowed and disallowed regions, respectively (Supplementary Table S3).

Residue conservation was reflected by surface mapping of the aa sequences of several plant TIPs over the generated EgTIP1.1 model (Figure 3). An overall evaluation of the EgTIP models detected the dual NPA (Asn-Pro-Ala) motifs conserved in all structures (Figures 3A,C; Supplementary Figure S1). In contrast, the key residues of the selectivity filter (positions H2, H5, LE1, and LE2) displayed some variations between the modeled proteins (Table 1; Figures 3B,C; Supplementary Figure S1). In all EgTIP1 and EgTIP2, the tetrad of the ar/R selectivity filter was composed of HIAV and HIGR, respectively (Table 1; Supplementary Figure S1). Despite showing phylogenetic relationships with TIP5 isoforms (Figure 1), the selectivity filter of EgTIP2.3 (HFGR) resembles the one found overall in group 2 members (HIGR; Figure 3E). Remarkably, EgTIP2.3 lacks the I residue in the H5 position, but has instead an F residue that preserves the hydrophobicity. Interestingly, the closest homologue of EgTIP2.3, VvTIP5.2, also harbors a classical TIP2-like selectivity filter (HIGR). For the other EgTIPs, slight residue variations at specific positions were observed (Table 1; Figure 3D). According to previous reports, a V residue is typically found at the LE2 position of TIPs from group 1, while an R is found in TIPs from groups 2, 3, and 4, whereas a C residue is present in all group 5 TIPs (Wallace and Roberts, 2004; Gupta and Sankararamakrishnan, 2009). In this regard, EgTIP3.1 harbors an atypical L residue in this position, while an R residue in EgTIP2.3 replaces the usual C residue (Table 1; Figure 3D). In addition, compared to the other EgTIPs, EgTIP3.2 and EgTIP5.1 harbor changes in the H2 position (S and N, respectively, instead of H), while EgTIP5.1 and EgTIP2.3 display variations in H5 (V and F, respectively, instead of I) (Table 1; Figure 3D).

Further comparison of the EgTIP-deduced aa sequences with those of other plant TIPs revealed that TIP members of groups 1 and 2 have the most conserved selectivity filters, including the H residue located in loop C (LC) of the proposed extended ar/R filter (Kirscht et al., 2016) (Figures 3D,E). VvTIP1.3 from V. vinifera, that possesses an M residue instead of a V in the LE2 position, is the only exception. The members of group 3 have most conserved residues in the H5 and LE1 positions and small variations in the H2 and LE2 positions. On the other hand, group 4 members have differences mainly in the LC and H5 positions, while EgTIP5.1 have the most variable LC positions.

To investigate the spatiotemporal expression profiles of the identified EgTIP genes, we first assessed their relative expression over a panel of different Eucalyptus organs/tissues using RT-qPCR. As shown in Figure 4, all members of the EgTIP subfamily, except EgTIP5.1, were transcribed. Among the genes investigated, solely EgTIP1.1 had a relatively broad distribution among organs/tissues, being expressed at moderate levels in the roots, stems, flowers and flower buds. In contrast, the other EgTIP genes surveyed showed distinctive organ/tissue-enriched expression. Of the genes grouped in the so-called group 1, EgTIP1.3 was preferentially expressed in vascular cambium and stems, EgTIP1.4 in flowers and flower buds, and EgTIP1.2 in leaves and stems. Of group 2, EgTIP2.1, and EgTIP2.2 were both moderately expressed in flowers and flower buds, but recapitulating previous data (Rodrigues et al., 2013), EgTIP2.2 was also highly expressed in roots. Of note, EgTIP2.3 was barely detected in leaves and stems and was not included in the heatmap. The duplicated gene pair EgTIP3.1 and EgTIP3.2 shared a similar fruit-enriched expression, but EgTIP3.2 was also detected at moderate levels in leaves. Finally, EgTIP4.1 was prominent in the stems.

EgTIP spatial expression patterns were further examined employing publicly available RNA-Seq data (FPKM values) from six different organs/tissues of a hybrid clone of Eucalyptus (Mizrachi et al., 2010). Interestingly, this dataset included vegetative and vascular tissue samples (such as phloem, xylem, and shoot tips) that were not represented in the previously performed RT-qPCR assays. According to the digital expression profiles generated (Figure 5), a wide expression pattern was observed for EgTIP2.1 and EgTIP1.1, which were detected at moderate to high levels in all organs/tissues examined. A similar pattern was observed for EgTIP1.2, except for its very low expression level in the phloem. According to our RT-qPCR assays, this gene was predominantly expressed in leaves (Figure 4). On the other hand, EgTIP1.4 was detected at moderate levels in mature leaves and at lower levels in young leaves and in shoot tips, while EgTIP4.1 was prominent in vascular tissues (phloem and xylem). This is consistent with the preferential expression of EgTIP4.1 in the stems as indicated by RT-qPCR (Figure 4). Contrasting with EgTIP5.1 that was not expressed, EgTIP2.3 was observed at very low levels in the phloem. The other EgTIP genes, including those mainly expressed in reproductive organs as determined by RT-qPCR, were almost undetectable in the RNA-Seq data surveyed. In this regard, the absence of EgTIP3.2 expression was particularly intriguing, especially because its transcripts were readily detectable in leaves by RT-qPCR (Figure 4). Overall, these results reveal a distinctive spatial distribution of the EgTIP genes, indicating that certain isoforms may act in a relatively organ-specific manner.

Inspection of the 5′ upstream regions (1200 bp) of the EgTIP genes revealed the presence of several putative cis-regulatory elements mainly implicated in phytohormone and abiotic stress responses (Table 2). In this context, all EgTIP promoter regions, with the exception of EgTIP1.2, contained more than one ABA-responsive element (ABRE), involved in hormone responsiveness and in ABA-mediated abiotic stress signaling (Yamaguchi-Shinozaki and Shinozaki, 2005). The presence of ABRE in the promoters investigated suggests a possible role of ABA in the control of EgTIP expression as already reported for HvTIP1.2 and HvTIP3.1 in barley (Lee et al., 2015). Another well-represented regulatory element found in all the promoters investigated, except EgTIP2.3, was the TGACG-motif implicated in plant responses to methyl jasmonate (MeJA), a well-known primary signal in plant defense responses. Most EgTIP promoters (9 out of 11) also contained a MBS (MYB-binding site) motif generally recognized by MYB transcription factors implicated in drought and salt inducible gene expression. The circadian and SP1 elements known to play a role in circadian control and in light responses, respectively, were also detected in the majority of the EgTIP promoters. The circadian element is similarly over-represented in the promoter regions of aquaporin genes from sorghum (Reddy et al., 2015). Additional stress-related regulatory elements found in a small number of EgTIP promoters include TC-rich motifs, LTR, HSE and C-repeat/DRE (Table 2). Overall, these findings indicate that the EgTIP genes are regulated by phytohormones such as ABA and MeJA and differentially responsive to a set of abiotic stimuli. This prediction is supported by a previous study showing the modulation of EgTIP2.2 promoter activity by osmotic stress and ABA (Rodrigues et al., 2013). However, the functionality of such predicted elements requires further validation.

Taking into account the aforementioned promoter analysis, we decided to evaluate the stress-responsive expression of each EgTIP gene. For this, total RNA extracted from roots, stems and leaves of Eucalyptus plantlets treated with PEG 8000 for 6, 12, and 24 h was employed. Relative to the untreated control, the majority of the investigated EgTIPs, with the obvious exception of the fruit-enriched EgTIP3.1 and EgTIP3.2 genes (not shown), showed increased expression upon PEG-induced osmotic stress at one or more time points (Figure 6). The timing and magnitude of the observed changes in gene expression were quite diverse among the investigated organs. In roots, most EgTIPs showed moderate fold-changes and reached peak expression levels mainly at 12 h of stress imposition (Figure 6A). In contrast, the timing of maximal expression in the stems varied between the investigated genes and was most evident at 6 and 12 h of PEG treatment (Figure 6B). EgTIP1.2 showed maximum stress-induced fold-change in this organ at 12 h of PEG treatment. Compared to roots and stems, EgTIP expression was increased in leaves. In this organ, significant fold-changes were observed at 6 h, and more substantially, at 24 h of stress imposition (Figure 6C). At this time point, EgTIP1.1, EgTIP1.2, EgTIP1.3, and EgTIP2.3 appeared as the most responsive genes. Interestingly, despite being barely detectable in leaves and stems, EgTIP2.3 was moderately and highly induced by osmotic stress in the investigated organs. In addition, the expression of EgPIP2, which was used as an inducible control (Rodrigues et al., 2013), was monitored in parallel and shown to be induced by PEG in all organs examined.

Based on our results, distinct organ/tissue expression profiles were observed for the EgTIP genes identified. Among them, EgTIP1.1 was the most widely expressed, indicating that it might play a more generalized role throughout the plant. This is in line with previous reports showing that TIP1.1 is widely distributed and the most abundantly expressed TIP isoform (Alexandersson et al., 2005; Knipfer et al., 2011; Regon et al., 2014). Conversely, the other EgTIP genes showed enriched-organ/tissue expression, a feature that was particularly evident in wood-related and reproductive organs/tissues. The only exception was EgTIP5.1, which appears to be nonfunctional since no transcripts or expression sequence tag evidences could be found for this gene. Interestingly, none of the identified EgTIP genes was found to be exclusively expressed in a specific organ/tissue.

Despite some partial overlaps, the observed divergences in their spatial expression patterns indicate that EgTIPs may have undergone functional diversification. This fact also suggests that certain isoforms might have acquired specific organ/tissue roles. Evidence from poplar trees, however, indicate that, in spite of their elevated tissue-specificity, a high level of functional redundancy exists amongst PtTIP members (Cohen et al., 2013). In this respect, the existence of functional redundancy between TIP genes was postulated in previous studies employing insertional Arabidopsis mutants (Schüssler et al., 2008). In general, current evidences indicate that TIP isoforms act redundantly to regulate water and solute transport. In contrast with this, Gattolin et al. (2009) demonstrated the existence of tissue- and cell-specific patterns of expression of some AtTIP genes in Arabidopsis roots. Similarly, a specific and developmentally regulated leaf expression pattern of some HvTIP subfamily members was observed in barley (Besse et al., 2011). These results highlight the existence of possible specialization of TIP isoforms, but further investigations are required to uncover possible functional specificity of the identified EgTIP genes.

In trees, the modulation of TIP expression in response to abiotic stress has been documented (Cohen et al., 2013; Rodrigues et al., 2013; Martins et al., 2015; Regon et al., 2014), and the major evidences point to a role for these proteins in stress adaptation. It is noteworthy that most of the EgTIP genes analyzed here displayed osmotic stress-responsive transcriptional profiles in three different organs. Although the timing and extent of their expression varied, these results indicate that the Eucalyptus responses to osmotic stress require the coordinated transcriptional regulation of most EgTIPs. The existence of cis-regulatory elements related to hormone and abiotic stress response in the EgTIP promoter regions lend further support to this idea. Consistently, we have previously shown the induction and repression of the EgTIP2.2 promoter activity by osmotic stress and ABA, respectively (Rodrigues et al., 2013). However, despite these reported evidences of transcriptional regulation, the possible contribution of the different EgTIP isoforms in the control of water flow under stress conditions requires further investigation.

The selectivity filter of aquaporins comprises four residues located at the non-cytosolic portion of the pore. These residues are implicated in water translocation and in the selection of other substrate molecules such as glycerol, ammonia, boric acid, CO₂, H₂O₂, and urea (reviewed in Maurel et al., 2015). Previous data have indicated that TIP subfamily members show the highest level of sequence divergence in the ar/R region among known plant aquaporins (Wallace and Roberts, 2004; Azad et al., 2012). Despite this, some residues are well conserved, particularly in the LE1 position, in which an A or G is present, and in the H₂ position, in which an H residue has been associated with water specificity (Sui et al., 2001; Kirscht et al., 2016). Remarkably, this conserved H residue is found in most identified EgTIPs, with the exception of EgTIP3.2 and EgTIP5.1.

Concerning LE2, an R residue known to contribute to the mechanism of proton exclusion is usually present at this position (Wu et al., 2009; Kirscht et al., 2016). In TIPs, however, this residue can be replaced by polar residues such as C or nonpolar residues such as V, a feature that is observed in some EgTIPs. Intriguingly, a nonpolar L residue in found in this position in EgTIP3.1. According to Azad et al. (2012), the presence of a hydrophobic residue in this position in TIP1 homologs renders the proteins multifunctional aqua-glyceroporins. Thus, changes in aa charge and side chains lead to modifications in pore size and hydrophobicity that have direct consequences in TIP transport selectivity.

A still unresolved question concerns EgTIP2.3. Due to its clustering with known TIP5 isoforms, EgTIP2.3 could be also assigned as a member of group 5. However, EgTIP2.3 lack the tetrad NVGC commonly found in the ar/R filter of most TIP5 isoforms. Instead, it harbors a selectivity filter that resembles the one present in TIP2 isoforms, an attribute that justifies its nomination as a TIP2 in Phytozome. The closest homologue of EgTIP2.3, VvTIP5.2, also harbors a canonic TIP2-like ar/R filter (HIGR). Intriguingly, such an ambiguous classification was also observed for certain NOD26-like intrinsic proteins from different plant species (Zou et al., 2015b). Overall, our results suggest that both EgTIP2.3 and VvTIP5.2 represent unusual isoforms sharing phylogenetic relationships with TIP5s but transport selectivity of TIP2s.

Several studies have demonstrated that additional residues may contribute to TIP preference for the transport of other substrates than water (Froger et al., 1998; Hove and Bhave, 2011). A recent study that elucidated the crystal structure of AtTIP2.1 proposed the expansion of the selectivity filter with the inclusion of a flexible H residue in loop C (Kirscht et al., 2016), that has been implicated in ammonia transport. In this respect, EgTIP2.1, EgTIP2.2, and EgTIP4.1 have the same H residue in loop C as AtTIP2.1. Comparatively, several TIP2 orthologs from other plant species also have an AtTIP2.1-like extended selectivity filter (Zou et al., 2015a,b, 2016), thus suggesting that these proteins may be able to transport ammonia. Moreover, in view of the predictions for non-aqua transport profiles proposed by Azad et al. (2016), both EgTIP2 and EgTIP4.1 as well as all EgTIP1 fulfill the requirements to act as H₂O₂ transporters.

Previous reports have also included the analysis of nine specificity-determining positions (SDPs) to predict non-aqua substrates (Hove and Bhave, 2011) as well as the Froger’s positions (P1-P5) to discriminate TIP-mediated glycerol transport (Froger et al., 1998). All these positions were evaluated in the EgTIPs investigated in order to address possible transport specificity. Although the Froger’s positions remained highly conserved in EgTIPs compared to TIPs found in other plant species (Table 1), they do not fit the requirements for glycerol transport. Noteworthy, as reported earlier for TIPs from Ricinus communis and Jatropha curcas (Zou et al., 2015b, 2016), the SDP analysis indicated the presence of residues typically associated with urea transport in all EgTIPs (Supplementary Table S4). In particular, EgTIP1.1, EgTIP4.1, and EgTIP5.1 have classical urea-type SDPs (Supplementary Table S4). Overall, these results are valuable for future studies on the structure–function relationships of the identified EgTIPs.