To identify the putative PA transporter family genes in barley, we used the known sequences of PUT and BAT genes of Arabidopsis and rice for BLASTp search (E value, 10⁻¹⁰) against the barley genome database hosted at Ensembl Plants (https://plants.ensembl.org/index.html). The protein sequences of barley were analysed using Hidden Markov Model (HMM) analysis for the presence of a typical HMM domain representing the corresponding protein class in a HMMER search (https://www.ebi.ac.uk/Tools/hmmer/) (Finn et al., 2011; Upadhyay and Mattoo, 2018). For accuracy, these sequences were also cross verified with InterProScan (Finn et al., 2017), and the candidates containing any of the typical domains were considered as proteins of interest.

Characteristics of identified proteins such as isoelectric point (pI), amino acid sequence length (aa), molecular weight (MW) and grand average of hydropathicity (GRAVY) were obtained using tools available at the ExPASY bioinformatics resource portal (https://www.expasy.org). The sub-cellular localizations were predicted by Plant-mPLoc webtool (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/). For gene structure analysis, genomic DNA and coding DNA sequences (CDS) corresponding to each identified gene were analysed. The GFF3/GTF annotation file containing the locations of interesting genes in genome and their structural information was extracted from the Ensembl Plants database (https://plants.ensembl.org/index.html). The MEME tool from the MEME suite 5.3.3 (http://meme-suite.org/tools/meme) was used to identify ten statistically significant motifs in the protein sequences based on “zero or one occurrence per sequence (zoops)” (Bailey et al., 2009). Hidden Markov Model analysis was done with HMMER database by the Inter-pro scan program hosted at webtool (http://www.ebi.ac.uk/interpro/). The gene structure, conserved motifs and domains were visualized using the TBtools software (Chen et al., 2018).

The protein sequences of reported PUTs and BATs in various plant species were retrieved and a phylogenetic tree was constructed using the Neighbor-Joining method with Poisson correction and 1000 bootstrap values using the MEGA-11 program (Tamura et al., 2021). The tree was visualized using the iTOLv6 program (https://itol.embl.de/). The evolutionary distances were computed using the p-distance method and are in the units of the number of amino acid differences per site (Tamura et al., 2021). All ambiguous positions were removed for each sequence pair (pairwise deletion option). Evolutionary analyses were conducted in MEGA-11 (Tamura et al., 2021). To explore the synteny relationships of the orthologous PA transporter genes among barley and other species, genome data and the gene annotation files of A. thaliana, O. sativa, B. distachyon, and S. bicolor were downloaded from the Phytozome database (https://phytozome-next.jgi.doe.gov/). The syntenic analysis graphs were constructed by using the Dual Synteny Plotter function in TBtools software.

To analyse the promoter regions of PA transporter genes, 1.5 kb upstream sequence from the translation start sites was retrieved from Ensembl Plants database. Transcription factor binding sites (TFbs) and cis-acting regulatory elements (CREs) were identified within the promoter regions to investigate their role in the regulation of polyamine transporter gene expression. The PlantPAN 3.0 software (http://plantpan.itps.ncku.edu.tw/promoter.php) was used to identify TFbs in the promoter regions of the PA transporter genes. The multiple promoter analysis program (http://plantpan2.itps.ncku.edu.tw/gene_group.php?#multipromoters) was used to identify the common TFbs in different promoter regions. The analysis of CREs was carried out using the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) and visualized with the TBtools software. The prediction of CpG/CpNpG islands and tandem repeats (TRs) in promoter regions was done by PlantPAN 3.0 web server (http://plantpan.itps.ncku.edu.tw/index.html). The coding sequences of barley polyamine transporter genes were analyzed by psRNATarget server (https://www.zhaolab.org/psRNATarget/) for miRNA target sites prediction (Dai et al., 2018). The protein–protein interaction analysis of PA transporters was carried out on the STRING web tool (https://string-db.org/) (Szklarczyk et al., 2021), and the network was generated using Cytoscape-3.9.0 software.

The prediction of secondary structure of HvBAT and HvPUT proteins was done using the SOPMA web server (Park and Igarashi, 2013). The three-dimensional (3D) structures of HvBAT/HvPUT proteins were predicted by using Phyre2 server (Seiler and Raul, 2005) and further evaluated with Ramachandran plot, ANOLEA (Atomic Non-Local Environment Assessment) and ProSA analyses. Additionally, the PROCHECK test was used to inspect the 3D structure of PA transporter proteins in the SAVES server (https://saves.mbi.ucla.edu/). The CLICK server was utilized to compare the protein models through the RMSD value calculation based on α-carbon superposition (Moinard et al., 2005). The active sites of HvBAT/HvPUT proteins were predicted through CASTp 3.0 server (Nakanishi and Cleveland, 2021). The post-translational modifications of HvPAT proteins were anticipated using MusiteDeep webserver (https://www.musite.net/) with default parameters. The 3D structures of ligands (PAs) viz., putrescine (C₄H₁₂N₂; PubChem ID: 1045), spermidine (C₇H₁₉N₃; PubChem ID: 1102) and spermine (C₁₀H₂₆N₄; PubChem ID: 1103) were downloded from PubChem database (https://pubchem.ncbi.nlm.nih.gov). The Autodock 4.2 and Autodock Vina softwares were utilized to prepare the receptor proteins and ligands, and docking simulations, respectively (Wang and Casero, 2006). Subsequently, HvBAT/HvPUT-PAs interactions were analysed with PyMOL (https://pymol.org) and visualised with LigPlot⁺ v.2.2.4 software (Mehta et al., 2002).

Initially, the identification of PA transporters included nine putative HvPUTs and seven HvBATs in barley genome; however two HvPUT and one HvBAT were excluded from further analyses due to their small/partial available sequences. Three out of seven HvPUTs (HvPUT2, HvPUT3 and HvPUT6) are located on chromosome 5, while chromosomes 3, 4 and 6 had one HvPUT gene each ( Table 1 ). All the HvBATs were located on chromosome 3 except HvBAT6, which was present on chromosome 2. The PA transporter proteins ranged in length from 328 (HvPUT5) to 527 (HvBAT6) amino acids, and their molecular weight varied from 34.85 kDa (HvPUT5) to 57.82 kDa (HvBAT6) ( Table 1 ). The theoretical pI ranged from 5.06 (HvPUT1) to 9.22 (HvBAT2), with an average of 8.04. Considering the stability of the proteins, three members (HvPUT1, HvPUT2 and HvPUT3) had the instability index of more than 40 while the rest of them were inferred to be stable with an instability index of less than 40. Moreover, all the PA transporter proteins had GRAVY values greater than 0, suggesting hydrophobic nature. The majority of HvPUT proteins had a high percentage of aliphatic amino acids, with an average aliphatic index of 104.87, ranging from 98.23 (HvPUT5) to 111.43 (HvPUT3). Furthermore, all the transporter proteins were predicted to be localized in the cell membrane and contained 12 transmembrane domains (TMDs), except HvPUT5 with seven TMDs ( Supplementary Figure S1 , Supplementary Table S1 ).

Upon further characterization of the PA transporters, we identified the distribution of ten conserved motifs in the studied proteins. The conserved motifs arranged in the order of motifs: 8, 2, 5, 1, 7, and 10 were shared by all HvPUTs except HvPUT5 which did not contain motifs 7 and 10 at the C-terminus ( Figure 1A ). Moreover, all the HvBAT proteins consisted of conserved motifs arranged in the order of 4, 2, 8, 6, 1, 10, 3, 5, and 9 ( Figure 1A ). All the PA transporter proteins consisted of the typical conserved domains for AA_permease_2 (PF13520) and AA_permease (PF00324) ( Figure 1B ). In addition, the PA transporters’ structural diversity was evaluated, and the distribution of exons and introns is shown in Figure 1C . Among HvPATs, most of the HvPUTs had one intron except HvPUT3 which consisted of three introns in the 5’UTR, while HvBATs consisted of 4-7 introns.

To establish the evolutionary and phylogenetic relationships of PA transporters, the protein sequences of identified PUTs and BATs from various plant species ( Supplementary Table S2 ) were utilized including: Arabidopsis (Arabidopsis thaliana), stiff brome (Brachypodium distachyon), rice (Oryza sativa), poplar (Populus trichocarpa), maize (Zea mays), wine (Vitis vinifera), saffron crocus (Crocus sativus), tomato (Solanum lycopersicum), and melon (Cucumis melo). The obtained result indicated that the PUT family divided mainly into three groups (Group I, II, and III) ( Figure 2 ). Group I included 12 proteins (including HvPUT1and HvPUT2), four proteins were classified into group II (including HvPUT2), and seven proteins were classified into group III which included HvPUT4-7 ( Figure 2 ). Similarly, BATs were also categorized in three groups; HvBAT1 and HvBAT2 were included into group IV, HvBAT3-5 were categorised in group V and group IV contained one HvBAT6. Overall, the HvPUTs and HvBATs showed close association with their homologs in other plant species.

To further understand the evolutionary relationship of PA transporters, collinearity analyses were conducted between barley, Arabidopsis, stiff brome, rice and sorghum (Sorghum bicolor) plants. A total of four PA transporters (30.76%) were identified with collinearity relationship in other monocot plant species ( Supplementary Figure S2 ; Supplementary Table S3 ). The four orthologous gene pairs (HvPUT1, HvPUT2, HvBAT3 and HvBAT6) were identified from H. vulgare-S. bicolor, H. vulgare-B. distachyon, and H. vulgare-O. sativa. Interestingly, no orthologous gene pair was detected from H. vulgare-A. thaliana.

The putative promoter regions (1.5 kb upstream sequences) of the PA transporter genes were analysed for the presence of cis-acting regulatory elements (CREs) and transcription factor binding sites (TFbs). A total of 1159 putative CREs were identified in the promoter regions of PA transporter genes. The identified CREs were further divided into four groups:160 elements were related to hormone response, 237 elements associated with stress response, 133 elements related to light response, and 629 elements associated with growth and development ( Figure 3 ; Supplementary Table S4 ). CREs for growth and development consisted of 19 elements, among which CAAT-box and TATA-box were the most abundant (288 and 251, respectively), and commonly shared by all of the transporter genes. Moreover, 21 different types of elements were found in the light responsive group, mainly: Box 4, G-box, A-box and CCGTCC elements. Comparatively, 15 types of elements were related to hormone response, mainly ABRE involved in abscisic acid (ABA) response, CGTCA-motif and TGACG-motif in jasmonic acid (MeJA) response, TGA-element and AuxRR-core in auxin response, and ERE in ethylene response. Additionally, stress responsive elements (19 types) comprised STRE (stress-response element), LTR (low temperature response), MBS for drought inducibility, WUN-motif for wound response, and ARE and GC-motif for anaerobic induction.

The putative promoter sequences of seven HvPUT and six HvBAT genes were used to find common TFbs types, namely TCP, WRKY, bHLH, NAC, BES1, bZIP, MYB, GATA and AP2/ERF. We identified 14,991 and 12,490 TFbs for HvPUT and for HvBAT, respectively ( Supplementary Table S5 ). In HvPUT and HvBAT gene promoters ( Figure 4 , Supplementary Table S5 ) the TFbs with the highest representation was the MYB-binding site (1151 and 1015, respectively) ( Figure 4 ), while the highest number of MYB TFbs was found in HvPUT2 (292), and HvBAT2 (249). The second most common TFbs was AP2/ERF-binding site, 988 for HvPUT (the highest number was in HvPUT3 with 366), and 875 for HvBAT (HvBAT3 with 258). The other common TFbs for HvPUT was TCP- (981) and bZIP-binding site (861), the highest numbers of these TFbs were present in HvPUT1 (196), and in HvPUT4 (201) promoter regions, respectively. While in HvBAT promoter regions the other common TFbs was bZIP- (787) and TCP-binding site (784), the highest numbers of these TFbs were present in HvBAT2 (175), and in HvBAT4 (183) promoter regions, respectively.

Furthermore, the presence of CpG/CpNpG islands and tandem repeats (TRs) were also analysed in the promoter regions of PA transporter genes. The CpG/CpNpG islands were identified in the 1.5 kb upstream promoter region of four HvPUT genes, namely HvPUT1, HvPUT3, HvPUT4 and HvPUT7, and of three HvBAT genes: HvBAT2, HvBAT3 and HvBAT4 ( Table 2 ). TR analysis showed that two HvPUTs (HvPUT6 and HvPUT7) contained TRs of 22 and 13 nucleotides in length, respectively, belonging to minisatellites while HvBAT4 has one TR of two nucleotides long and belong to microsatellites ( Table 3 ).

A total of six barley miRNAs (hvu-miR) comprising target sites in the CDS regions of seven PA transporter genes were identified ( Figure 5 ; Supplementary Table S6 ). Among HvPUTs, The HvPUT7 consisted of two target sites for hvu-miR6189 and hvu-miR6196, while HvPUT2, HvPUT3 and HvPUT6 had target sites for hvu-miR6181, hvu-miR1120 and hvu-miR6196, respectively. Moreover, HvBAT1 and HvBAT2 had target sites for hvu-miR6185, whilst HvBAT5 contained target site for hvu-miR6191. Almost all of the identified miRNA-targeted genes, except for HvPUT7, were predicted to be silenced by cleavage inhibition. To represent the target accessibility, energy required to unpair the secondary structure around target site (UPE) was also calculated. The results indicated that the UPE varied from 8.657 (hvu-miR6189) to 24.211 (hvu-miR6196) ( Supplementary Table S6 ).

In silico analyses of the protein–protein interactions (PPI) demonstrated that only five HvPUT proteins (HvPUT1, 2, 3, 4 and 7) out of all 13 analyzed (seven HvPUTs and six HvBATs) show protein–protein interactions ( Figure 6 ). All the HvPUTs (1, 2, 3, 4 and 7) showed interactions with proteins A0A287F0C4 (Proteasome subunit beta), A0A287GHA7 (26S proteasome non-ATPase regulatory subunit 1 homolog), A0A287J5L6 (PCI domain-containing protein), A0A287RUQ0 (Proteasome subunit beta), A0A287RXT6 (Proteasome subunit beta type), A0A287V2I9 (AAA domain-containing protein), and F2D894 (Proteasome subunit alpha type). Moreover, HvPUT1-3 showed interactions with proteins A0A287DZ26 (DEK_C domain-containing protein), A0A287L3B8 (ABC transporter domain-containing protein) and A0A287NSR3 (DEK_C domain-containing protein). Furthermore, A0A287PCD4 (VWFA domain-containing protein) protein had interactions with HvPUT1, 2, 3, and 7 ( Supplementary Table S7 ).

Based on our secondary structure analysis of PA transporters, the predominant structure of the HvPUT proteins was alpha helix (42.99-53.91%), followed by random coil (28.46-38.72%), extended strand (13.03-14.75%), and beta turn (3.81-5.94%) ( Supplementary Table S8 ). Moreover, the alpha helix in HvBAT proteins ranged from 43.91% to 47.43%, the random coil varied from 30.27% to 33.72%, and the extended strand ranged from 18.03% to 20.50% ( Supplementary Table S8 ). Surprisingly, beta turn was only present in HvPUT proteins. The three-dimensional structure of PA transporter proteins was created using the best available template; the confidence in 3D modeling of HvPUT and HvBAT proteins was 100%, while coverage with the best template varied from 79% (HvPUT5) to 88% (HvPUT6) ( Figure 7 ; Supplementary Table S9 ).

Ramachandran plot analysis of 3D structures showed that most of the PA transporters had >85%, <16.7% and <1% of residues distribution in most favoured, favoured and disallowed regions, respectively ( Supplementary Figures S3 , S4 ). However, HvPUT5 had 76.2%, 22.4% and 1.4% of residues distribution in most favoured, favoured and disallowed regions, respectively. The ProSA Z-scores of barley PA transporters structures ranged from −3.42 (HvPUT2) to -1.0 (HvBAT2) and the ANOLEA Z-scores ranged from 7.56 (HvBAT6) to 12.22 (HvPUT5) ( Supplementary Table S10 )

We further carried out the analysis for prediction of post-translational modification (PTM) sites in PA transporters of barley. In HvBAT and HvPUT proteins, nine possible PTMs (acetylation, glycosylation, hydroxylation, SUMOylation, methylation, phosphorylation, palmitoylation, ubiquitination, and pyrrolidone carboxylic acid) were predicted ( Supplementary Figure S5 ; Supplementary Table S11 ). Most of the PTMs were randomly distributed in the proximity of the N-terminus of all proteins except of HvPUT2 which had PTMs mainly close to C-terminus. The most common PTMs included nine hydroxylation sites (HvPUT1) followed by six phosphorylation sites in HvPUT7 protein, five acetylation sites in HvPUT4 protein, and four glycosylation sites in HvPUTs (1,3,5,6, and 7) and HvBAT (1 and 3) proteins. The SUMOylation site was only predicted in HvPUT3 and HvPUT4, whereas one palmitoylation site was observed in HvBAT (1 and 5) proteins.

To explore the interactions of PA molecules in the binding site of predicted PA transporters, we used a molecular docking method. PA transporters showed the lowest mean binding energy (ΔG) for spermine (Spm) (ΔG: −4.8 kcal/mol; pki: 3.56 μmol) compared to spermidine (Spd) (ΔG: −4.5 kcal/mol; pki: 3.36 μmol) and putrescine (Put) (ΔG: −3.8 kcal/mol; pki: 2.79 μmol). Among all the HvPUTs, the HvPUT7 showed the lowest ΔG with Put (ΔG: −4.4 kcal/mol; pki: 3.22 μmol) and Spm (ΔG: −5.6 kcal/mol; pki: 4.11 μmol), while HvPUT6 had the lowest ΔG with Spd (ΔG: −4.9 kcal/mol; pki: 3.59 μmol). Moreover, HvBAT3 showed the lowest ΔG with Spd (ΔG: −4.9 kcal/mol; pki: 3.59 μmol) and Spm (ΔG: −5.2 kcal/mol; pki: 3.81 μmol) whereas, HvBAT5 showed the lowest ΔG with PUT (ΔG: −3.8 kcal/mol; pki: 2.78 μmol) ( Table 4 ; Figure 8 ).

In terms of HvPUTs-PA interactions, the number of formed hydrogen bonds varied from one in HvPUT1-Put, HvBAT1-Spd, HvBAT1-Spm and HvBAT2-Spm interactions to six in HvBAT3-Spd and HvPUT3-Spm interactions ( Supplementary Table S12 ). HvPUT7 with the lowest ΔG showed four hydrogen bonds and four hydrophobic interactions with putrescine, and three hydrogen bonds and ten hydrophobic interactions with spermine ( Figure 8 ). HvPUT6 had three hydrogen bonds and eight hydrophobic interactions with spermidine. In case of HvBATs, HvBAT5 with the lowest ΔG showed two hydrogen bonds and nine hydrophobic interactions with putrescine. Similarly, HvBAT3 with the lowest ΔG showed six hydrogen bonds and five hydrophobic interactions with spermidine, and showed four hydrogen bonds and 11 hydrophobic interactions with spermine ( Figure 8 ).

To study the possible involvement of specific PA transporters in different phases of development, we have performed gene expression analyses in barley plants. Our studies demonstrate that HvPUT1-4 show relatively stable and high expression level during plant growth and development stages namely: germination, seedling, tillering, stem elongation, booting, flowering, milk ( Figure 9A ). HvPUT5 was expressed only at seedling stage, while HvPUT6 was expressed in stem elongation-, flowering-, and milk stages. Expression of HvPUT7 was not detected in barley plants at any stage of development. Moreover, HvBAT1, 2 and 6 were stably expressed during almost whole plant development, while HvBAT3 and HvBAT5 were mainly expressed in seedling- and flowering stage, and HvBAT4 was detected at very low level exclusively in seedlings ( Figure 9A ). Unsurprisingly, all PA transporter genes show very low or undetectable expression in the last stage of barley plant development, namely dough stage.

The organ specific expression analysis ( Figure 9B ) revealed that among HvPUTs, only HvPUT3 is expressed in high level in almost all studied organs, while HvPUT2 and HvPUT4 showed high expression in most of the organs, except caryopsis and its elements (i.e embryo, endosperm, aleurone layer and scutellum) as well as in the shoot apex (HvPUT4). The other HvPUTs showed very low (HvPUT1, 5 and 6) or non-detectable (HvPUT7) expression in the studied plant organs. Among HvBATs, HvBAT1 and 2 showed high expression in the most of analyzed organs; the exception were caryopsis and its elements, as well as microspore (HvBAT1) and leaf (HvBAT2). Interestingly, HvBAT3, 5 and 6 were expressed mainly in roots and seedlings (although the expression of HvBAT3 was on a very low level), while HvBAT4 was not detected in any studied organs ( Figure 9B ).

Additionally, the expression of PA transporter genes in barley was analysed under various abiotic stress conditions ( Figure 10 ). HvPUT2 was clearly upregulated in wounding stress, and downregulated during cold stress, while the other HvPUT genes did not show any significant expression changes in response to abiotic stresses. Moreover, wounding stress resulted in clear upregulation of HvBAT1 and 2 genes as well as in a slight increase of HvBAT5 and 6 expression. Simulated drought stress conditions led to decreases in HvBAT5 and 6 transcript level in barley roots, and was accompanied by upregulation of the HvBAT4 gene. Moreover, temperature stresses seem to effect HvBATs expression in leaves, leading to upregulation of HvBAT1,2 and 6 in cold (mainly in prolonged cold-treatment experiments) followed by downregulation of HvBAT3 and 6 in heat-drought combined stress. Salt treatment showed minor effect on PA transporter genes expression level ( Figure 10 ).

The RNA- sequencing data developed in our laboratory allowed us to analyze changes in PA transporter gene expression during leaf senescence in barley plants. Clear downregulation of HvPUT5 expression was observed during the dark-induced leaf senescence (DILS) experiment, while the other HvPUTs did not show any significant changes in their expression ( Figure 11 ). Moreover, the expression of HvBAT3, HvBAT6 and HvBAT5 was upregulated in both DILS- and developmental leaf senescence (DLS) samples. HvBAT2 had slightly higher expression specifically in DLS and control conditions. The expression of HvBAT1, was downregulated during DILS, however, no change was observed in DLS.

In this work, seven PUTs (HvPUT1-7) and six BATs (HvBAT1-6) genes were identified as putative PA transporters in the barley genome. The average theoretical pI value, and the GRAVY values suggest that the majority of PA transporter genes encode alkaline-, hydrophobic- and more likely membranous proteins, that is additionally supported by the subcellular localization predictions as all the PA transporters were localized to the cell membrane. The motif analysis showed that all the PA transporters shared highly conserved corresponding sequences, indicating their similar functionalities. HvPUT5 did not contain motifs 7 and 10 at the C-terminal region. This may be due to partial protein sequence or HvPUT5 may have undergone motif deletion. The protein motif deletion or rearrangement facilitates functional protein diversification, contributing to flexible protein network evolution (Saito et al., 2007). The sequence alignment, protein domains, and motif analysis indicate that the PA transporter genes have shared the same evolutionary pattern and have similar function. All the PA transporters proteins had the signature HMM AA_permease_2 (PF13520) and AA_permease (PF00324) confirming the amino acid transporter activity. The PA transporters’ structural diversity was evaluated by assessing the distribution of exons and introns. The variation in the number of introns is anticipated because the number and length of introns in genes vary depending on organism and gene structure, and these differences may be related to intron function (Grabherr et al., 2011). Moreover, the gain and loss of introns can change the structure of genes and play a vital role in the evolution of gene families (Xu et al., 2012). HvPUT7 and HvBAT4 had the largest intron, it can be assumed that these genes served as a precursor from which the other members originated by undergoing an ‘intron loss’ event through the time of evolution (Gilbert et al., 1997). Also, HvPUT5 has a large intron but a very short second exon, which may confirm the conjecture made earlier that this may be due to partial protein sequence. We note that HvPUT3 had three introns in the 5’UTR which depicts its uniqueness among all the HvPUT members. We hypothesize that this gene structure might be involved in the regulation of gene expression or RNA stability (Shi et al., 2020).

Our results divided the PA transporters in different groups and the barley PA transporters were clustered together with other plants transporters implying that they may serve similar functions to their plant homologs. The barley genome had a different number of homologs for PA transporters as compared to Arabidopsis implying that the barley PA transporter gene subfamilies may have undergone gene expansion and/or gene loss during the course of its evolutionary history.

To further investigate the evolutionary relationships of PA transporter genes, we performed genome-to-genome synteny analysis between barley and four representative plant species. Comparatively, the PA transporter genes in Arabidopsis did not show any orthologous correlation with barley. However, the monocotyledons showed a considerable synteny. Therefore, we speculated that the syntenic correlations between PA transporters might be connected to the species’ evolutionary divergence. In particular, four PA transporters were identified to be syntenic across all the monocotyledons tested, indicating that these orthologous pairs are conserved and may have existed before the species divergence. We note that the intersections of syntenic PA transporters among different species may be useful for undertaking future relevant gene evolution studies.

We investigated the potential regulatory mechanism that controls the expression of PA transporter genes in barley by looking into cis-acting regulatory elements (CREs), transcription factor binding sites (TFBs) and CpG/CpNpG islands in the promoter regions. Here, we identified a total of 1159 putative CREs taking part in multiple biological processes and found a high number of CAAT-box and TATA-box elements in the promoter regions of PA transporters. Additionally, a large quantity of light-, hormone-, and abiotic stress-responsive elements in the promoter regions were identified which imply that the expression PA transporters might be controlled by a number of different elements which act through several pathways. We further identified the TFbs and classified as common TFbs (TCP, WRKY, bHLH, NAC, BES1, bZIP, MYB, GATA and AP2/ERF). These results further support the CREs analysis that expression of PA transporters might be controlled by a number of different elements. For instance in our study, the MYB-binding site was the most abundant TFbs in all barley PA transporter gene promoters. The highest number of MYB TFbs was found in HvPUT2 gene promoter, therefore, it may be strong candidate for MYB TFbs-mediated barley development control.

The presence of CpG/CpNpG islands in promoter regions represents critical sites for DNA methylation, which resulted in gene silencing (Cao and Jacobsen, 2002; Zemach and Grafi, 2007). The CpG/CpNpG islands were identified in the promoter region of four HvPUTs and three HvBATs, implying that expression of these genes might be controlled by DNA methylation. Overall, the study of CREs, TFbs and CpG/CpNpG islands suggested that promoters of PA transporter genes may play important roles in modulating the gene expression during developmental processes and/or stress tolerance. Further, we analysed the tandem repeats (TRs) in the promoter regions. In the promoter sequences of the genes we studied, TRs were found only in HvPUT6, HvPUT7 and HvBAT4. This indicates that in these genes there is a higher probability of mutation accumulating during replication (known as polymerase slippage) (Ahmar and Gruszka, 2022). The presence of TRs in these promoters may be useful for mutational analysis, and it may have a role in controlling gene expression of PA transporters.

Many reports have shown that miRNAs regulate responses of barley to different stress conditions (Kantar et al., 2010; Hackenberg et al., 2012; Sunkar et al., 2012; Wu et al., 2018). In this study, six barley miRNAs with target sites in CDS of seven PA transporters were identified. These miRNAs might be involved in various stress responses in barley. For instance, upregulation of hvu-miR6196 occurs during salt adaptation of the autopolyploid Hordeum bulbosum (Liu and Sun, 2017) and during barley exposure to an excess of boron (Unver and Tombuloglu, 2020). In this study, hvu-miR6196 had target sites in two transporters; HvPUT2 and HvPUT7 implying that their expression might be regulated by miRNA during the various stress conditions. Thus, exploring the role of miRNAs in PA transporter gene functions in response to various stresses would be of great interest.

Protein–protein interactions (PPI) represent an important aspect of plant systems biology which enables uncovering the unknown functions of proteins at the molecular level and can provide insight into complex cellular networks. In this study, the regulatory PPI network for the barley PA transporter proteins indicated considerable interactive networks among HvPUTs and several other proteins. This result implies that PUTs take part in protein complexes which have essential roles in regulatory processes, cellular functions and signaling cascades. For instance, all the HvPUTs showed PPI with proteasome subunits (alpha/beta). Interestingly, HvBAT protein did not show any PPI, suggesting that HvBATs are not involved in the regulatory protein complexes which have essential roles in physiological processes. Alternatively, the possibility exists that the database for PPI network of BATs has not been fully established making this type of interpretation more challenging.

In this study, we carried out homology modeling of the PA transporter proteins, and the 3D structures of all the proteins showed the high accuracy of the structures predicted. ProSA Z-scores analysis revealed the structures were typically plotted within the range of scores for native proteins of similar size from X-ray crystallography and NMR sources, suggesting no significant deviation from the native structures. Furthermore, we predicted the post-translational modifications (PTMs) in the protein sequences of barley PA transporters. Protein PTMs significantly increase proteome diversity, increase functionality, and allow for rapid responses, all at a low cost to the cell. The commonly found PTM sites were glycosylation and phosphorylation and were present in all the PA transporters. Our results suggested that barley PA transporters might play various roles in development and growth, and various stress conditions through PTMs. Interestingly, some PTM sites were found on specific transporters such as SUMOylation sites were only predicted in HvPUTs (3 and 4), palmitoylation sites were observed in HvBATs (1 and 5), suggesting their specific roles in plant physiology.

Probing the PAs binding pockets of barley PA transporter proteins through molecular docking facilitated an in silico understanding of the mechanisms and interactions involved in HvPUT/HvBAT-mediated PAs transport. Weak intermolecular interactions such as hydrogen bonding and hydrophobic interactions are key players in stabilizing energetically-favored ligands, in an open-conformational environment of protein structures (Patil et al., 2010). The numbers of hydrogen bonds in the PA transporter-PA complex were less than the number of hydrophobic interactions in all protein-ligand complex. Our results showed that the most suitable ligand-protein interaction included HvPUT7-Put/Spm, HvBAT3-Spm/Spd, HvBAT5-Put, and HvPUT6-Spd. However, further experimental analysis should be done to verify these predictions. Overall, these results provided that the interacting residues are functionally conserved and have important considerations for protein activity while designing PA transporters modified by genetic manipulation in barley.

The organ-specific expression analysis showed that mRNA of three out of the seven identified barley PUT genes, namely HvPUT2, HvPUT3, and HvPUT4, were clearly detectible in a majority of plant organs. The observation that three listed above may be the main transporters responsible for PA uptake in plants is additionally supported by the relatively stable expression pattern of these three genes during plant growth and development. Interestingly, HvPUT1 and HvPUT6 genes characterized by low expression in different barley organs, show transient upregulation during flowering and early seed development suggesting their specific functions in these processes. Among all PA transporters, HvPUT7 was rarely expressed in all the analysed organs and developmental stages. The archived down-regulation in some gene expressions is fundamental for keeping up with the gene duplicates and ancestral functions (Qian et al., 2010). Thus, the archived down-regulation of HvPUT7 expression in barley might be regarded as critical for maintaining their biological processes and protecting them from any misfortune during cell crisis. Our study showed that HvBAT3 and HvBAT5 were mainly expressed in seedling and flowering stage, suggesting their role in these stages. HvBAT1, HvBAT2 and HvPUT4 showed high expression in most of the organs except endosperm, shoot apex, embryo, caryopsis, aleurone layer and scutellum.

Changes in the PA pools in plants under stress conditions are well documented and have been reviewed (Wang et al., 2003; Wang and Casero, 2006; Nakanishi and Cleveland, 2021; Pál et al., 2021). Several PA biosynthesis genes were highly upregulated by cold, salt and/or drought stress and, as a result, cellular PA concentrations were increased in Arabidopsis under these stress conditions (Urano et al., 2003; Hummel et al., 2004; Urano et al., 2004; Vergnolle et al., 2005; Alcázar et al., 2010). The question remains whether the expression of PA transporters may be also altered in response to stresses in order to fine tune the optimal PA concentration in the cell or to improve intra/extra cellular transport rates to maintain homeostasis. In our analysis, HvPUT2 was upregulated in wounding stress, while downregulated during cold stress, suggesting its role in these stresses. HvBAT3 was upregulated in cold stress and HvBAT4 had high expression under simulated drought stress. HvBAT5 showed higher expression in salt stress, and HvBAT6 was downregulated in drought, heat;drought and simulated drought stresses.

Furthermore, in this study we analysed the expression of PA transporters during leaf senescence (DILS and DLS) in barley. Our results showed that HvPUT4 was slightly and HvPUT5 was strongly downregulated in DILS. The expression of HvPUT1 and HvPUT2 has not changed during DILS. The expression of HvBAT3, HvBAT6 and HvBAT5 was upregulated in both, DILS and DLS. PUTs are mainly responsible for PA uptake and BATs may act as exporters of PAs, as mentioned before. Thus, we speculated that in senescing cells the uptake of PAs was halted and they mainly relied on the biosynthesis of PAs, which is also supported by findings that PA biosynthetic genes are upregulated during leaf senescence in barley (Sobieszczuk-Nowicka et al., 2016; Tanwar et al., 2022). PA transporters might play important roles in the regulation the leaf senescence by regulating cellular PA homeostasis. The upregulation of HvBATs suggest that PAs may have possible roles in N remobilization to other cells/organs.

The molecular processes underlying leaf senescence are strongly conserved between plant species, suggesting that senescence has evolved as a selectable trait in plants. The phenomenon of senescence is often portrayed as a paradox, as this trait promotes the death of the individuals. This view is too simplistic, as plants are not slated to die before they undergo successful reproduction. Senescence is essential for sustaining the phenotypic plasticity of growth, and it represents an important evolutionary trait that enables plants to adapt to the environment. The impact of senescence on crop yield and quality, and its potential use in breeding more environmentally-resilient plants, are becoming increasingly important. An important aspect of that is the adequate uptake and remobilization of N that increases the plant’s N usage efficiency, thereby reducing the requirement for fertilizers (Schippers et al., 2015; Paluch-Lubawa et al., 2021). Since cells cannot store reduced N as NH₃ or NH₄ ⁺, PAs are important N storage compounds in plants (Wuddineh et al., 2018).