Genome-wide identification of Fagus sylvatica aquaporins and their comparative spring and summer expression profiles

A total of 45 aquaporins was identified in Fagus sylvatica, 35 of which were differentially expressed in spring and summer in the leaves, phloem and xylem of 11-year-old trees. European beech (Fagus sylvatica) has been widely studied in terms of its water relations and local adaptation. However, to date, the underlying basis conferring adaptation to differences in water availability are unknown. Therefore, we examined the expression of aquaporins in trees of four different beech provenances representing the southern and northern range margins, as well as core populations, grown in a common garden. We sampled their xylem, phloem and leaf tissue, when leaves had fully expanded, and in late summer. A total of 45 aquaporin isoforms were identified in the beech genome, of which 35 were detected across all sampled tissues. In our phylogenetic analysis, beech aquaporins clustered into the five subfamilies found in other woody species. Members of the plasma membrane intrinsic protein subfamily generally displayed the highest levels of expression, followed by tonoplast intrinsic proteins. Isoforms of the remaining subfamilies, Noduline-26-like intrinsic proteins, small basic intrinsic proteins and uncharacterised intrinsic proteins, were expressed at very low to moderate levels. The expression of most isoforms was stable or declined from spring to summer. Leaves followed a different expression profile from that of vascular tissues, whereas both phloem and xylem were found to express the same FsMIPs. Tissue-specific aquaporin expression was very similar amongst the four beech provenances, indicating that there is no inherent difference in the capability of these provenances to regulate aquaporin activity. The general decrease in FsMIP expression toward the end of the growing period indicates that aquaporins are involved in tree water relations and growth.


Introduction
European beech (Fagus sylvatica) is one of the most widely distributed forest tree species, spanning multiple climatic zones from the Mediterranean to Scandinavia. It is shadetolerant and generally susceptible to drought (Robson et al. Communicated by Resco de Dios. 2012; Gömöry et al. 2015), but it is genetically highly diverse due to its biogeographical history (Bontemps et al. 2016). This genetic diversity forms the basis for the high phenotypic plasticity observed in European beech (von Wühlisch 2008;Robson et al. 2012;Sánchez-Gómez et al. 2013). Common-garden trials at its southern range margin have found that local adaptations also confer stress tolerance, allowing beeches from less favourable southern sites to withstand moderate droughts (Robson et al. 2012;Sánchez-Gómez et al. 2013;Gömöry et al. 2015). These differences in the sensitivity to drought, and thus water use, of beech populations can have wide implications for future forest management in a changing climate with summer drought predicted to increase in frequency and severity (IPCC 2022).
In trees, two systems function in parallel to transport water and nutrients to and from the leaves. In the xylem, water flow is driven by transpiration from the roots and through the leaves, according to the cohesion-tension theory (Dixon and Joly 1894). The driving force for the movement of solutes in the phloem is more debated, but it is generally accepted that the turgor gradient generated by the loading of assimilates in the leaves' phloem lowers the osmotic potential and draws water in, while unloading of assimilates in the sink tissues raises osmotic potential and draws water out to the surrounding tissue (Münch 1930). Phloem and xylem are located in close proximity, and hydraulically connected by radial cells in the wood tissue, which leads to exchanges of water between the two compartments (Hölttä et al. 2006). In the leaves, stomatal movements modulate both tension in the xylem via transpiration and phloem transport via sugar loading capacity by regulating CO 2 assimilation (Cowan and Farquhar 1977). Therefore, it would be physiologically beneficial to regulate water transport jointly in both tissues, in accordance with leaf physiology. Combined models of xylem and phloem water flows have proved useful to elucidate their hydraulic link (Hölttä et al. 2006), but the specific mechanisms of regulation, and genotypic variability that would produce diverse stress responses, have yet to be discovered. Now that the F. sylvatica genome has been sequenced (Mishra et al. 2018(Mishra et al. , 2021, it has become more feasible to examine the molecular basis for the hydraulic link between phloem and xylem as well as the role of specific genes in relation to water use strategies. Thus, determining the genetic basis for differences in fitness among populations under specific stresses has become a realistic goal. One strong candidate for conferring increased resistance to drought stress are aquaporin genes, as aquaporins have been firmly established to play an important part in various processes related to plant water balance (Zimmermann et al. 2004;Secchi and Zwieniecki 2014;Secchi et al. 2017). In trees, aquaporins have been implicated in diameter changes of the stem throughout the day (Steppe et al. 2012), embolism refilling (Secchi and Zwieniecki 2014;Secchi et al. 2017), long-distance water translocation (Zimmermann et al. 2004) and adjusting root pressure (Wegner 2014).
Aquaporins belong to the super family of major intrinsic proteins (MIP) and are divided into five subfamilies, as shown in Table 1. Aquaporin isoforms in tree species, are more diverse than isoforms of the more extensively studied model and crop species which lack the entire subfamily of X intrinsic proteins (XIP, see references of Table 1). Genome duplication events, as well as gain or loss of exons and introns, have caused the diversity of aquaporins to drastically increase in some species, such as poplar (Populus sp.) and cotton (Gossypium sp.), while purging events resulted in the loss of XIPs in others, such as the monocots (Gupta and Sankararamakrishnan 2009;Sonah et al. 2017;Jiang et al. 2020;Venisse et al. 2021). In addition, aquaporins, like any other gene, are subject to natural selection pressures, resulting in differences in their numbers across species as well as differences in their relative abundance across the five subfamilies (Table 1). Despite their name and the earliest research on them, aquaporins are not entirely water-specific, but may conduct a range of different substrates from CO 2 and H 2 O 2 to urea and ammonia (Kojima et al. 2006;Bertl and Kaldenhoff 2007;Hooijmaijers et al. 2012;Uehlein et al. 2012). The substrate specificity of an aquaporin depends on two main features, which ultimately determine the potential roles of an isoform in plant water relations. First, the narrowest parts of the channel, namely, the signature sequence NPA and the ar/R constriction prevent large molecules from passing through while at the same time excluding protons and removing single water molecules from their hydration shell, thus ensuring a high rate of water movement across the membrane (Ilan et al. 2004;Wallace and Roberts 2004;Tornroth-Horsefield et al. 2006;Horner et al. 2015;Kitchen et al. 2019). Second, aquaporin permeability to substrates other than water is less dependent on pore size itself, than on nine amino acid residues lining the pore (Zou et al. 2015a;Kitchen et al. 2019). These residues are, therefore, called substrate determining positions (SDP). SDPs can easily be identified through the alignment of amino acid sequences and allow the substrate permeability of newly discovered, yet uncharacterised aquaporins, to be predicted. Substrate prediction, therefore, facilitates the evaluation of whether a newly discovered aquaporin isoform expressed in a certain tissue could play a role in regulating water homeostasis, and whether this may be accomplished through the direct impact on water movements or indirectly through the accumulation of solutes.
Despite the broad range of permeabilities displayed by the many aquaporin isoforms, they all share a common feature: the capability to conduct water across membranes. On the level of the whole plant, the role of aquaporins in water balance has been the focus of considerable attention, but with only moderate success in obtaining clear identification of individual isoforms even under controlled conditions (Voicu et al. 2009;Ben Baaziz et al. 2012;Secchi and Zwieniecki 2014;Kromdijk et al. 2020;Israel et al. 2021;Wang et al. 2021b). This gap in our knowledge has made it very difficult to interpret ecophysiological data, especially for trees that are often studied under much less controlled conditions than model plants. Though their precise roles remain speculative, aquaporins are nevertheless believed to be involved in a variety of processes related to tree water balance.
To evaluate the role of aquaporins in regulating water use in different provenances of F. sylvatica, we identified the complete set of aquaporins they possess, as well as their gene expression levels in the xylem, phloem, and leaves for four provenances spanning the entire species distribution range. We compared four distinct provenances in a common garden which had been well-characterised in terms of their water relations and performance (Wang et al. 2021a). Our aim was to test for differences in aquaporin gene expression between spring, when xylem vessels were still differentiating but leaves had fully expanded, and late summer, when vascular tissues were fully differentiated. To our knowledge, this is the first comparative examination of expression changes between xylem and phloem tissues in the bark of trees.

Plant material and sample collection
Beech seedlings from four provenances representing the south-westerly range margin (Montejo, Spain), core distribution range (Eichelberg, Germany), a high elevation environment (Rindelpholz, Germany) as well as the northern range margin (Blaviksliarna, Sweden) were planted in a common garden experiment in Helsinki (N 60 13′, E 25 01′) in autumn of 2010. Details on the climate at the sites of origin and treatment of beechnuts are given in Wang et al. (2021a). In spring and late summer 2021, four 11-year-old trees per provenance were sampled during the growing season (14.-16.06.2021 and 30.-31.08.2021).
Samples for RNA sequencing were taken in the morning between 9 and 11 a.m. when photosynthetic capacity and transpiration, and thus water movement through the tree, were close to their daily maximum. Thus, we aimed to detect transcription of aquaporin genes when they were likely to be most active in water movement (Steppe et al. 2012). Healthy and fully expanded leaves were chosen to be sampled from a height between 2 and 3 m above ground, from a branch receiving direct sunlight for at least some part of the morning hours. The sampled leaves were immediately frozen in liquid nitrogen. Separate phloem and xylem samples were collected using a scalpel to peel off a strip of bark from each tree trunk (dimensions ca. 2 × 8 cm) approximately 50 cm above ground. Thin layers of xylem from the trunk, and phloem from the inside of the bark piece, were scraped off following the procedure described in Sakr et al. (2003) and quickly dried using blotting paper, after which they were immediately frozen in liquid nitrogen.
The frozen samples were ground in Eppendorf tubes inserted into Teflon blocks pre-cooled to − 75 °C using a shaker (Retsch MM301, Verder Scientific/Retsch, Düsseldorf, Germany) for 3 min at 30 cycles per second. On average, 50-80 mg of the ground powder was used for RNA extraction. RNA was extracted using the GeneJET Plant RNA Purification Mini Kit (ThermoFisher Scientific, MA, USA) according to the manufacturer's instructions with the exception that the Plant RNA Lysis Solution was supplemented with β-mercaptoethanol instead of DTT, as well as an additional 100 µl Plant RNA Isolation Aid (Ambion Inc., TX, USA) to remove polysaccharides and polyphenols from potentially lignified tissues. The lysis solution with the added sample was homogenised for 20-30 s and centrifuged at 12,000×g for 5 min. The supernatant was collected and used to continue normally with the RNA extraction according to the manufacturer's instructions. The quality and concentration of the extracted RNA was determined using a ND-1000 Spectrophotometer (Ther-moFisher Scientific) and samples with low yield (< 10 ng/ µl) or low purity (A260/280 < 1.7 and/or A260/230 < 1.8) were discarded, resulting in n = 3-4 for each provenance. All RNA samples were treated with DNase1 to remove potential contamination with genomic DNA. RNA sequencing was performed by the Functional Genomics Unit (FuGu) at the University of Helsinki. In total 33 952 genes were identified across all samples.

Identification of aquaporin genes and bioinformatics analysis
A BLASTp search using the Arabidopsis MIPs identified a total of 54 potential FsMIPs, including pseudogenes, genes transcribed into identical proteins and splice variants, in the F. sylvatica genome (beechgenome.net, Supplemental Table S1) (Altschul et al. 1997;Mishra et al. 2021). After removing pseudogenes and splice variants, the remaining 45 full-length protein sequences were aligned with ClustalW (Thompson et al. 1994) in the MEGA11 software (Tamura et al. 2021) together with the 35 AtMIPs. A phylogenetic tree was constructed using the neighbourjoining method and branch reliability tested with 1000 bootstrap replicates.

Gene expression analysis
Read count tables from the RNA sequencing were analysed in R with the package DESeq2 (Love et al. 2014). p values indicating statistically significant up-/down-regulation of aquaporin isoforms were Benjamini-Hochberg adjusted to control for false discovery rate (Benjamini and Hochberg 1995). A principle component analysis (PCA) was performed including all 33,952 detected genes, to remove potential sample outliers which could have resulted from very low gene counts due to insufficient RNA quality. After excluding outliers, n was 3-4 for each sample time-tissueby-provenance combination. Gene expression analysis was continued using the 35 identified FsMIPs detected in our samples, to determine potential differences between provenances. No significant differences were found among provenances, and therefore, the four provenances were pooled in the subsequent analysis resulting in n = 12-15 for each sample time-by-tissue combination. A second PCA was performed in R with the package FactoMineR  on the normalised gene counts for all 35 FsMIPs detected in our samples, and a biplot was generated from this analysis using the package factoextra (v1.07). The heat map was generated with the R package gplots (v3.1.1). All artwork was finalised in CorelDRAW Graphics Suite 2017.

Results
We identified 45 full-length non-redundant aquaporin genes in the F. sylvatica genome, among which we found 12 PIPs, 12 TIPs, 12 NIPs, 3 SIPs, and 6 belonging to the XIP family that is not found in Arabidopsis Venkatesh et al. 2015). In line with previous reports on other species, PIPs, SIPs and XIPs can be further divided into two sub-groups, whereas TIPs and NIPs show a much higher diversity. The phylogenetic tree including the FsMIPs and AtMIPs clearly illustrates the division into the four subfamilies common to both species and that the XIPs form a separate fifth cluster (Fig. 1). The basic protein properties of all FsMIPs, their residues forming the ar/R selectivity filter and the predicted range of non-water substrates is summarised in Table 2, together with their proposed nomenclature. The identified proteins vary in length from 181 to 429 amino acids and their genomic distribution is uneven across all chromosomes, which is consistent with earlier reports (references listed in Table 1). Some of the MIPs described in Table 2 have fewer predicted transmembrane units than the six that would be expected for putative aquaporins and these MIPs also generally lack one of the NPA signature sequences, but were nevertheless expressed in our samples.
Gene structure analysis (Fig. 2) revealed that FsPIPs are very homogenous with typically 4 exons and 3 introns, FsTIPs are more variable with 2-4 exons and 1-3 introns and FsNIPs display the highest diversity with 3-7 exons and 3-6 introns. With an average of 2 each, FsSIPs and FsXIPs had the lowest numbers of introns and exons. This genetic structure is very similar to that found among aquaporins in other plant species (Chaumont et al. 2001;Johanson et al. 2001;Sakurai et al. 2005;Gupta and Sankararamakrishnan 2009;Park et al. 2010;Zou et al. 2015b;Liu et al. 2019;Wei et al. 2019;Venisse et al. 2021). Furthermore, conserved gene motifs as identified in MEME, are also shown in Fig. 2 and exemplify the highly conserved structure of aquaporins in general, and within each subfamily in particular.
In addition to predicting the number of transmembrane units using the TMHMM-2.0 online tool, we used the structurally well-characterised SoPIP2;1 as a template for the identification of transmembrane units as well as residues of the ar/R selectivity filter ( Figure S1 and Table 2). FsPIPs, belonging to the same subfamily as SoPIP2;1, share a high sequence identity and are thus very similar in structure with 6 transmembrane units as well as identical ar/R constrictions. Alignment with SoPIP2;1 also clearly shows that FsPIP1;5 possesses 6 transmembrane units despite the fact that the online tool predicted only 4 such units. Similar results were obtained for FsTIPs, FsNIPs and FsSIPs with the exception that in this case those isoforms predicted to have fewer transmembrane units were lacking the corresponding amino acid sequences. However, as these aquaporins were nevertheless detected in our samples, and must,   therefore, be functional genes/proteins, we assumed incomplete genome sequencing to be the reason behind the truncated nature of the predicted proteins. FsXIPs were found to be similar in structure to the other subfamilies, but the predicted transmembrane units TM2 and TM4 were longer and all FsXIPs contained the three conserved cysteins common to XIPs of other plant species (Zou et al. 2015a).
Substrate specificity determining positions (SDPs), from previously identified aquaporins with permeabilities to substrates other than water, were aligned with FsMIPs to predict their potential range of permeabilities. The SDPs identified for all 45 FsMIPs, as well as their predicted substrates, are summarised in Tables 2 and S2. Multiple FsMIPs were predicted to be permeable to boric acid, hydrogen peroxide as well as urea including all FsPIPs. FsPIPs, with the exception of FsPIP2;6 were also predicted to conduct CO 2 . Interestingly, as the only non-PIP, also FsXIP2;1 may be permeable to CO 2 . Aquaporins of the FsNIP3 subfamily were predicted to conduct silicic acid, while only one FsMIP, FsTIP2;1, may conduct ammonia. Based on the identification of SDPs, none of the FsSIPs appeared to be permeable to substrates other than water.
Of the 45 identified FsMIPs, 35 were detected in our samples (Table 2, Figs. 2, 3). With the exception of FsNIP1;1 in leaves, expression of all FsNIPs was very low. FsNIP1;3 and FsNIP7;1 were leaf-specific, while FsNIP3;1 was only expressed in bark tissues. Most of the highly expressed aquaporins belonged to the FsPIPs. FsPIP1;3 and FsPIP2;2 were the most highly expressed isoforms in xylem and phloem, while FsPIP2;3 was the most abundantly expressed isoform in leaves. Expression of FsPIP1;4, FsPIP1;5 and FsPIP2;5 appeared to be leaf-specific, though their abundance in leaves remained low throughout the trial period and extremely low levels of expression were also found in phloem. Of the three FsSIPs, only FsSIP1;2 was expressed at moderate levels in all sampled tissues, but all FsSIPs were detected at both sampling times and in all tissues. FsTIP1;1 and FsTIP4;4, were the most abundantly expressed TIPs with high levels of expression in xylem and phloem, but very low levels in the leaves. The overall expression of FsTIPs was very low in leaves. Furthermore, FsTIP3;1 was found to be leaf-specific, FsTIP4;1 xylem-specific and FsTIP5;1 was not detected in leaves. Only one FsXIP, FsXIP2;1, was detected in our samples: it was expressed at very low levels and restricted to the xylem and phloem.
In terms of the variation in gene expression between spring and late summer, almost all FsMIPs found in our samples were detected at both timepoints in all tissues investigated (Fig. 3). Expression levels of almost all isoforms decreased in late summer compared to the spring (Figs. 3,  4). In the leaves, the expression of few genes changed significantly, whereas we found a much higher variability in phloem and xylem (Table 3). Most of the statistically  connecting loop E *Divergence due to incomplete sequencing significant changes in expression occurred for lowly expressed genes, with highly expressed genes maintaining their expression levels between spring and late summer. The log2-fold change for all 18 significantly up-/down-regulated MIPs is given in Table 3. Six of these genes were more highly expressed in at least one tissue type in late summer, 13 were expressed less, while the expression of all remaining FsMIPs not listed in Table 3 did not significantly differ in their expression between spring and summer in any of the sampled tissues.
Expression levels of FsMIPs differed most strongly in leaves between spring and summer, with only FsTIPs being expressed at consistent levels (Fig. 4). Phloem and xylem differed considerably in their spring and summer expression profiles from that of leaves, and experienced less pronounced changes in expression levels between the seasons (Fig. 4). Although their profiles were similar, phloem and xylem nevertheless differed in the genes that were increased/ decreased from spring to summer. In phloem, only FsPIP2;4 increased its summer expression compared to spring, while FsPIP1;2, FsTIP1;3, FsTIP1;4, FsTIP2;1 and FsXIP2;1 decreased. In xylem, an increase was observed in FsPIP1;1 and FsSIP1;2, with the decreased expression levels found in FsPIP2;6, FsPIP2;7, FsTIP2;2, FsTIP4;2 and FsTIP4;4. The expression of two genes, FsPIP1;3 and FsTIP1;4, changed in a similar manner in both tissues, with FsPIP1;3 having increased in summer and FsTIP1;4 decreased.
A principal component analysis (PCA) of all 35 FsMIPs detected in our samples revealed a close similarity in gene expression between phloem and xylem tissues, whereas leaf samples formed a separate cluster (Fig. 5). PC1 describes variation between phloem/xylem and leaf, and thus variation between tissues, and accounts for 23.3% of the total variation as apportioned by the PCA. PC2 accounts for 10.5% of total variation and it segregates each tissue into distinct clusters according to sampling time, with the spring samples located on the positive side of the axis and summer samples on the negative side. Overall, FsMIPs are more strongly expressed in spring than late summer as can be seen from their contribution to the biplot variables (Fig. 5). In leaves, a total of eight FsMIPs (FsNIP1;1, FsPIP1 ;4,FsPIP1;5,FsPIP2;3,FsTIP1;2, FsTIP4;1; FsSIP1;1 and FsSIP2;1) contributed Shown is the mean of the gene counts separated according to tissue type ± SE. The SE was too small for most genes to be vis-ible in the graph and has, therefore, been omitted. Asterisks indicate statistically significant up-/down-regulation in late summer compared to spring Fig. 4 Heatmap of the mean FsMIP expression in the sampled tissues. Letters L, P and X indicate leave, phloem and xylem samples, respectively, whereas numbers 1 and 2 indicate spring and summer samples, respectively. The row z-score is based on the mean expression of a particular FsMIP across the tissues, providing an estimate on the level of gene expression for each sample, and is displayed for each sample as a deviation of the black solid line from the dashed mid-line (mean). Expression levels below the average row z-score are shown in red and expression levels higher than the row z-score in white. The dendrogram is based on the column means of the z-score strongly to the expression profile in the spring, all of which were diminished in late summer (Figs. 3, 4, Table 3). Only FsSIP1;2 experienced a significant increase in expression in leaves in late summer. The expression profile in phloem/ xylem in the spring was dominated by five FsTIPs (1;1, 1;3, 1;4, 2;1 and 4;4) and one FsPIP (2;6, Fig. 5). This colour-coded as indicated by the legend. Larger symbols in the same colour as the individual measurements indicated the mean for the spring (open symbols) and late summer (closed symbols) sampling located on the positive and negative side of the PC2 axis, respectively. Black symbols denote the mean for the respective tissue overall. The contribution of the individual genes to the tissueand season-specific variation is given by the distance from the origin as indicated by the length and transparency of the connecting lines, with a long and dark line indicating a large contribution TIP-dominated expression shifted toward a more PIP-dominated expression profile in late summer, namely, through the increased expression of FsPIP1;1 and FsPIP1;3. In overall, however, phloem and xylem experienced a reduction in MIP expression in late summer, with FsSIP1;2 being the only gene to display a significantly increased contribution (Figs. 3, 4, Table 3).

Discussion
We identified 45 non-redundant aquaporin isoforms in the F. sylvatica genome, which, in light of the 33 isoforms discovered in silver birch (Betula pendula) and the 55 isoforms in poplar (Populus trichocarpa), is well within the expected range (Gupta and Sankararamakrishnan 2009;Venisse et al. 2021). Likewise, the distribution of the identified isoforms across the five subfamilies and their overall genomic structure in terms of the exon-intron arrangement is as we would expect, when compared to other woody species (Wei et al. 2019;Jiang et al. 2020;Venisse et al. 2021).
The 45 isoforms that we identified include several incomplete gene sequences (Tables 2, S1), which is not unusual among recently sequenced woody species (Park et al. 2010;Wei et al. 2019). The predicted proteins of these incomplete sequences were either uncharacteristically small, missing one of the NPA-signature sequences, transmembrane units or residues of the ar/R constriction. However, the RNA sequencing detected these genes in our samples, which lead us to the conclusion that they must be fully functional aquaporins, as they would not be expressed otherwise. The partial amino acid sequences of these genes, furthermore, account for the low bootstrap values in parts of our phylogenetic analysis.
The NPA sequences were very well-conserved within the FsPIP and FsTIP subfamilies, as well as for most FsNIPs, with FsSIPs and FsXIPs showing the highest degree of divergence, which is in agreement with reports on other species (Johanson et al. 2001;Gupta and Sankararamakrishnan 2009;Park et al. 2010;Wei et al. 2019). Likewise, the ar/R selectivity filter of FsPIPs was found to be very well-conserved within the subfamily as well as across species, whereas FsTIPs and FsNIPs showed more variability within the subfamily, but to a similar extent as in other plant species. FsSIPs and FsXIPs were found to be the most variable subfamilies in terms of their ar/R residues (Gupta and Sankararamakrishnan 2009).
Originally identified to be water-specific channel proteins in plants (Maurel et al. 1993), the range of substrates able to cross membranes via aquaporins has subsequently greatly expanded. Among these, for example, plasma membrane intrinsic proteins (PIP) have been found to also be permeable to CO 2 , H 2 O 2 and Na + (Heckwolf et al. 2011;Hooijmaijers et al. 2012;Uehlein et al. 2012;Byrt et al. 2017;Groszmann et al. 2021); tonoplast intrinsic proteins (TIP) conduct ammonia (NH 3 ) and urea (CH 4 N 2 O) in addition to water (Kojima et al. 2006;Bertl and Kaldenhoff 2007), while  show the broadest range of substrate specificities from boric and silicic acid to ammonia, glycerol, urea and lactic acid (Wallace et al. 2002;Takano et al. 2006;Choi and Roberts 2007;Zou et al. 2015a). XIPs, as the most recently discovered subfamily have not been studied as intensively, but are known to conduct boric acid, H 2 O 2 and urea (Zou et al. 2015a). As such, small and basic intrinsic proteins (SIP) are the only aquaporin subfamily for which non-water substrates have not been reported and they may thus be truly water specific.
Substrate prediction using SDPs identified a range of potential non-water substrates with the FsPIPs displaying the highest diversity including boric acid, CO 2 , H 2 O 2 and urea, which is in line with previous studies (Bienert et al. 2014;Zou et al. 2015b). Due to the high degree of conservation within PIPs, these substrates were predicted to permeate almost all identified isoforms. CO 2 permeability is characteristic of and specific to PIPs (Wei et al. 2019). However, according to our prediction, FsXIP2;1 is also a potential conductor of CO 2 . FsNIPs were the second subfamily that we found to include a large range of predicted non-water substrates including boric acid, H 2 O 2 , silicic acid and urea. Silicic acid permeability is a common feature of NIP2s across many species and can easily be identified by the highly conserved ar/R selectivity filter of this subgroup (Bansal and Sankararamakrishnan 2007;Chiba et al. 2009;Park et al. 2010;Zou et al. 2015a;Wei et al. 2019). The most commonly predicted non-water substrate for FsTIPs was urea, but FsTIP2;1 was predicted to be also permeable to ammonia, which is a typical feature of TIP2s (Zou et al. 2015b;Wei et al. 2019). None of the FsSIPs were predicted to have any non-water substrate permeability, which is consistent with the lack of existing reports of such permeability (Zou et al. 2015a), and in light of their localisation to the endoplasmic reticulum  perhaps not surprising. However, the ar/R constriction of FsSIPs and SIPs in general is variable, which would suggest SIPs to be permeable to substrates other than water (Venisse et al. 2021). Interestingly, FsXIPs also had remarkably few potential non-water substrates, which is surprising, because XIPs have been reported to conduct urea and boric acid (Wei et al. 2019) as well as glycerol and H 2 O 2 (Deshmukh et al. 2013).
There were in total 35 FsMIPs expressed in the collected samples, most of which were found in all three tissue types in spring as well as late summer, confirming the ubiquitous nature of aquaporins (Fraysse et al. 2005). In agreement with reports that PIPs are abundant in the plasma membrane, FsPIPs were among the most highly expressed isoforms (Robinson et al. 1996;Johansson et al. 2000;Monneuse et al. 2011;Byrt et al. 2017) in all three tissues Kirch et al. 2000;Secchi et al. 2017). This PIPdominated expression was particularly apparent in leaves in spring, with high expression levels of FsPIP1;3, FsPIP2;2 and FsPIP2;3. This might suggest that these three isoforms play a role in transpiring and photosynthetically active leaves, as reported for PIPs in other species (Heckwolf et al. 2011;Ben Baaziz et al. 2012;Israel et al. 2021), and may also support their predicted CO 2 permeability (Table 2). It is noteworthy that PIP1s with PIP2s are co-expressed, not only in leaves, but phloem and xylem as well (Fig. 3). PIP1s require the interaction with PIP2s in the tetrameric assembly to be exported from the endoplasmic reticulum and localise to the plasma membrane (Fetter et al. 2004;Zelazny et al. 2007;Otto et al. 2010;Groszmann et al. 2021). However, PIP2s have also been reported to modulate their activity through interactions with each other . This is likely to occur in our samples as well, because we detected three highly expressed PIP2 isoforms in each sampled tissue and almost all PIP2 isoforms identified in this study were also present in the three tissues (Fig. 3).
The high expression levels of FsPIP1;3 and FsPIP2;2 in the vascular tissues suggest that these two isoforms are involved in water transport (Stanfield et al. 2017). This could be over long distances, or between phloem and xylem as well as the surrounding tissue, as the covariation of FsPIP1;3 and FsTIP1;4 between spring and summer may indicate. Water exchange between phloem and xylem is dynamic and occurs to balance out changes in xylem pressure caused by fluctuations in rates of transpiration (Patrick et al. 2001;Fraysse et al. 2005). Interestingly, the expression of FsPIP1;3 in both the xylem and phloem increased rather than decrease toward the end of the summer; likewise, FsPIP1;1 and FsPIP2;1. In late summer, at a time when overall rates of transpiration are expected to decrease, increased expression of these three PIP isoforms may support the hypothesis that PIPs in bark tissue principally facilitate radial water exchange between tissues, cycling of water between phloem and xylem, or embolism repair rather than long-distance water translocation (Sevanto et al. 2011;Steppe et al. 2012;Secchi and Zwieniecki 2014;Secchi et al. 2017).
Aquaporin expression in phloem as well as xylem was not as PIP-dominated as that in leaves, because of relatively high contributions of FsTIP1;1 and FsTIP4;4, especially in the phloem. However, irrespective of their localisation to the tonoplast, TIPs have been reported to be highly expressed in vascular tissue, to facilitate water and solute diffusion between tissues, and to regulate hydraulic conductance (Beebo et al. 2009;Pou et al. 2013). In this context, their high phloem/xylem expression in spring when cell expansion occurs is not surprising as TIP1s are involved in cell division and elongation (Ludevid et al. 1992;Chaumont et al. 1998). It may also point to FsTIPs largely determining the hydraulic conductance between phloem and xylem and, therefore, radial water exchange (Sevanto et al. 2011). Sevanto et al. (2011) suggested that the hydraulic coupling of phloem and xylem tissue is very tight, which would require a large radial hydraulic conductance to balance rapidly changing water fluxes. This would explain why we detected not only high levels of FsPIP expression in these two tissues but also high levels of FsTIP1;1 and FsTIP4;4.
Despite the similarities in their expression profiles, phloem differed from xylem in the genes that changed in expression in the summer (Fig. 4). Since our study is the first comparative examination of expression changes between xylem and phloem tissues in the bark of trees, supporting evidence in the literature is scarce, which forces us to some degree of speculation. Both tissues mature over the course of the summer and this is likely to involve the differential regulation of aquaporin genes. Many FsTIPs displayed lower expression levels in summer, when cell expansion and lignification has ceased, which points to FsTIPs being likely regulators of xylem and phloem development. Similar conclusions can be drawn for FsPIP2;6 and FsPIP2;7 in xylem, as their expression decreased significantly in the summer when rates of transpiration are expected to decrease. Both isoforms are predicted to conduct H 2 O 2 , which is required for lignification (Marjamaa et al. 2007). Therefore, their expression may be needed in newly developed xylem tissue during the spring. This would also explain lower expression levels during the summer, as the tissue matures. Phloem differed considerably from xylem in the strength and significance of their expression change in the summer. FsPIP1;1 and FsPIP1;3 increased in expression in summer in the phloem, a change seen in xylem as well, though not as pronounced. This change may have led to an increased radial water exchange between the tissues. By contrast, the expression of FsPIP1;2 was significantly lower in summer in phloem, but not in the xylem. FsPIP1;2 is, therefore, more likely to contribute to phloem cell expansion and maturation, than radial water exchange.
Overall, expression of FsMIPs was highest in the spring, when trees were actively growing and transpiring and their contribution decreased toward late summer. This shift in expression levels is reflected in the distinct clusters formed by the spring and summer samples in our principal component analysis (Fig. 5) and could at least in part be related to the lower temperatures in late August, since aquaporin expression and activity is temperature sensitive (Jang et al. 2004(Jang et al. , 2007Lee et al. 2012;Venisse et al. 2021). We did not find any differences in aquaporin expression between the beech provenances during either spring or summer. Therefore, there appears to be no intrinsic difference in the capability for these provenances to regulate their water balance through changes in aquaporin activity when grown in the climate of southern Finland. However, since beech provenances vary in their susceptibility to stresses, such as drought (Robson et al. 2012;Sánchez-Gómez et al. 2013;Gömöry et al. 2015), it is conceivable that genetically diverse provenances differ in their regulation of aquaporin activity when exposed to abiotic stress, as such differences in regulation have been reported for drought in cultivars of Citrus sinensis (Wei et al. 2019), Malus domestica (Liu et al. 2013) and Phaseolus vulgaris (Zupin et al. 2017).

Conclusions
We identified 45 Fagus sylvatica aquaporins; a similar diversity to that found in other studied tree species, and they cluster into the same five subfamilies identified previously. Aquaporins of the PIP subfamily were among the most abundantly expressed in leaves as well as vascular tissues and one of their roles may be to facilitate radial water exchange. We also found the overall aquaporin expression to decrease differentially between spring and late summer, supporting the hypothesis that they play a major role in regulating tree water relations especially at times of high transpirational demand and high rates of photosynthesis. The predicted substrate permeabilities of FsMIPs point toward potential physiological roles of the various isoforms, thus knowledge of their differential expression can be useful in selecting more stress-resistant tree populations for forestry in a changing climate.
Author contribution statement DI, YS and TMR conceived the ideas for the manuscript. DI and TMR designed the experiments. DI designed the field methodology. DI collected and analysed the data with significant contributions from MD. DI and MD wrote the manuscript with editorial contributions from all co-authors. TMR supervised the work with contributions from JJZ and YS.
Funding Open Access funding provided by University of Helsinki including Helsinki University Central Hospital. This research was funded by Academy of Finland decisions #324555, #323843 and #304519, and partially supported by two awards from Viikki Plant Science Centre (ViPS) and HiLIFE Grand Challenges Programme at the University of Helsinki granted to YS and TMR.

Data availability
The full RNA sequencing data set generated during the current study is available from the corresponding author on reasonable request. The raw RNA sequencing data containing the identified FsMIPs analysed in this study is available in Supplementary file S2.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.