The PSB strains used in this experiment are part of a collection reported in our previous study (Koczorski et al., 2022). The bacteria were isolated (autumn 2018) from the roots and rhizosphere soil of two willow species: ‘Tora’ (Svalöf-Weibull (SW species No. 910007, S. schwerinii × S. viminalis) and ‘Loden’ (SW 890129, S. dasyclados). These willow trees were grown at two test sites established in 2014 as part of the ECOLINK project in Uppsala (Sweden) and Rostock (Germany) (Hoeber et al., 2018; Koczorski et al., 2021b). The selection of PSB was performed using three selection media containing calcium triphosphate (National Botanical Research Institute’s phosphate growth medium (NBRIP) and Pikovskaya (PVK)) or calcium diphosphate (DCP) as phosphorus sources (Koczorski et al., 2021a). The bacteria that showed the ability to solubilize phosphorus (presence of halo zones) were classified into the PSB group and identified based on their 16S rRNA sequences (OP102593-OP102680) (Koczorski et al., 2022).

In the present study, a total of 64 bacterial strains from the PSB collection belonging to the species Pseudomonas, Bacillus, Erwinia, Serratia, Paenibacillus or Burkholderia (Koczorski et al., 2022) were screened. Solid NBRIP and DCP media were used to select the phosphate-solubilizing bacterial strains with the highest efficiency for phosphorus solubilization [the composition of the media is given in (Koczorski et al., 2021a)]. The bacteria were spot inoculated onto media (4 strains per Petri plate), and there were three replicates. The bacteria were cultured in the dark at 28°C, and the diameter of the halo zone was measured after 1, 5 and 10 days of culture. The results obtained were statistically analyzed; the 10 most efficient phosphate-solubilizing bacterial strains were selected and used for further analysis. The 10 selected bacteria did not include strains that were pathogenic to humans and plants, despite their high activity (data verified based on published literature).

The ten phosphate-solubilizing bacterial strains selected in the previous step were tested in liquid NBRIP media. NBRIP medium was chosen since plants cannot efficiently solubilize TCP on their own (Ticconi and Abel, 2004). In this part of the experiment, three different phosphorus concentrations were used, namely, 0 mM (NP), 1.0 mM (0.3129 g/L TCP) (LP) and 2.0 mM (0.6259 g/L TCP) (HP), and noninoculated media with the same variants were maintained as controls. The concentration of available phosphorus in the liquid NBRIP medium for all 10 strains was determined using the spectrophotometric molybdenum blue method (30 samples in total) (Nagul et al., 2015). Based on a statistical analysis, the 2 most effective PSB strains, Paenibacillus spp. and Pantoea agglomerans, were selected and used to inoculate the plants in the pot experiment.

Pots were divided into two compartments using membranes and filled with sterile sand (pH around 7). Loden (SW 890129, S. dasyclados) and Tora (SW 910007, S. schwerinii × S. viminalis) willow cuttings with two or more nodes were inserted into the pots such that there were two cuttings per pot (total 6 plants per variant). To minimize the effect of phosphorus sorption in the soil, quartz sand was used as a substrate. The pot experiment was designed to include three phosphorus concentrations and two willow species for inoculation treatments (control (Ctr), bacterium 1 (Pantoea agglomerans accession number: OP102615) and bacterium 2 (Paenibacillus spp. accession number: OP102612). Both bacteria were previously isolated from surface sterilized roots of Tora species. Prior to inoculation bacteria were cultivated on R2A medium (BD Difco, USA) for 2 days in 28°C then diluted to 0.7 OD600 in sterile water right before inoculation. The temperature in the room was set to 22°C with a photoperiod comprising 10 hours of light and 14 hours of darkness. After the cuttings acclimated for 7 days, with the exception of the control plants, the pots were inoculated with 10 ml of bacteria inoculum. The plants were watered thrice a week with a fertilizer solution with one of three different phosphorus concentrations (NP, LP, HP) for 1 month. We aimed to reach the final phosphorus concentrations (0.3129 g/L for LP and 0.6259 g/L for HP) in the first week of the experiment. After the plants were cultivated for 5 weeks, we harvested them and measured the length of their roots and shoots, along with measuring the fresh weights of their roots, shoots and leaves. In addition, the total phosphorus content in leaves at the end of the experiment was examined using soil fractionation method (Negassa et al., 2010) and via inductively coupled plasma–optical emission spectrometry (ICP–OES). The available phosphorus in soil substrate was measured using the spectrophotometric molybdenum blue method after sample extraction with 1% citric acid solution (Van Reeuwijk, 2002) for all the analyzed variants in the experiment (180 samples in total). A portion of the leaves obtained at harvest was flash frozen in liquid nitrogen and stored at -80°C.

The plant growth parameters (leaf weight, shoot weight, root weight, shoot length, root length) and phosphorus concentrations in leaves and soil were analyzed using STATISTICA v.13.3.721 (StatSoft, Poland). As the results obtained were not normally distributed, the Kruskal−Wallis test was used to determine significant differences. Dunn’s post hoc test was performed to determine differences between inoculation variants, different phosphorus concentrations and the willow genotypes used. The data were analyzed using 3-way ANOVA, with soil phosphorus concentration (3 levels), Salix species (2 levels) and bacterial inoculation (3 levels) included as factors, and their interactions were also evaluated. The response variables included all the growth parameters and phosphorus concentrations in the soil and plant leaves.

Total RNA was extracted from 100 mg of leaf tissue and crushed in sterile mortars filled with liquid nitrogen. RNA extraction was performed according to the protocol of Chomczynski and Sacchi (2006), with slight modifications. A total of 1000 µl of TRIzol reagent (Life Technologies, Poland) was added to the crushed leaf material, which was then incubated for 5 minutes at room temperature. Then, 300 µl of chloroform was added, after which the material was incubated for 15 minutes and then centrifuged at 12000 rpm for 15 minutes. A mixture of 0.8 M sodium chloride, 1.2 M sodium citrate and 200 µl of isopropanol was added to the upper layer of the supernatant, after which the mixture was centrifuged. The resulting pellet was washed twice with 75% ethanol. The RNA pellet was then resuspended in RNAse-free water (50 µl) and stored at -80°C. The RNA samples were outsourced to the company Novogene (Great Britain) for library preparation and transcriptome analysis.

For transcriptome analysis each variant of the experiment was represented by three samples. Messenger RNA was purified using poly-T oligo-attached magnetic beads. The first strand cDNA was synthesized using random hexamer primers, followed by the second strand cDNA synthesis using either dUTP for directional library or dTTP for non-directional library. The library was checked with Qubit and real-time PCR for quantification and bioanalyzer for size distribution detection. Quantified libraries will be pooled and sequenced on Illumina platforms, according to effective library concentration and data amount. The total RNA was sequenced using the Illumina platform in conjunction with P5 (P5-AATGATACGGCGACCACCGAGA (5’-3’) and P5’-TTACTATGCCGCTGGTGGCTCT (3’-5’) and P7 adapters (CGTATGCCGTCTTCTGCTTG-P7’ (5’-3’) and CATACGGCAGAAGACGAAC-P7 (3’-5’). After the removal of reads that were contaminated with adapter sequences, that included unknown nucleotides that constituted more than 10% of either read (N > 10%) and that were of low quality (base quality < 5) and constituted more than 50% of the sequence, the genome of Salix viminalis retrieved from the European Nucleotide Archive (ENA) was used as a reference for mapping using HISAT2 (Almeida et al., 2020). Assembly of the sequences was performed using StringTie software with transcripts of class code type ‘u’. Quantification was performed using featureCounts, and a differential analysis was performed using DESeq2 and edgeR with the following parameters: |log2(fold-change)| >= 1 &padj<= 0.05. Finally, enrichment analysis was performed using clusterProfiler with a padj< 0.05.

The solubilization indices obtained from the experiment (given in the range from 0 to 1.6) are presented in the form of a heatmap where for each experimental variant (growth medium: DCP and NBRIP, strain: 64 in total and time of cultivation: 1, 5 and 10 days), the values are indicated based on color differences (from blue to yellow) ( Figures 1A, B ). The values in Figure 1 are sorted according to the phosphorus solubilization efficiency measured on the NBRIP medium after 1 day of bacterial cultivation as the main criterion for the selection of microorganisms for further studies.

Of the 64 strains tested, 15 with the strongest ability to solubilize phosphorus on DCP and NBRIP solid media were selected. In this group, we observed a predominance of species of the genus Bacillus (7 out of 15: Bacillus megaterium, Bacillus aryabhattai, Bacillus spp. and Bacillus mycoides); Erwinia (2 out of 15: Erwinia spp.); and single strains of Rhizobium spp., Streptomyces spp., Pantoea agglomerans, Paenibacillus spp., Caballeronia glathei and Serratia spp. The strains with the highest activity on NBRIP media after 10 days of incubation belonged to Bacillus megaterium, Erwinia spp., Bacillus spp. and Bacillus aryabhattai. The second DCP medium showed a slightly different distribution of strains in terms of their activity. Among the 15 strains previously selected (based on the NBRIP medium), only four displayed high activity for phosphorus solubilization of DCP: Caballeronia glathei, Bacillus aryabhattai, Pantoea agglomerans and Paenibacillus spp. Notably, the strains with the highest activity on the DCP medium included Pseudomonas frederiksbergensis, Mitsuaria spp., Paenibacillus spp. and Pseudomonas spp. Among the group of strains, Bacillus megaterium, Bacillus aryabhattai, Erwinia spp. and Caballeronia glathei were excluded from further study, as they were previously reported to be potential plant pathogens (Abdel-Monaim et al., 2012; Narsing Rao et al., 2019, Dobritsa and Samadpour, 2019, Parcey et al., 2022). A more detailed quantitative determination of bacterial phosphorus solubilization efficiency in liquid media (NBRIP) with different phosphorus contents showed that Rhizobium spp., Pantoea agglomerans, and Paenibacillus spp. were most efficient at solubilizing phosphorus compounds ( Figure 2 ). Among the strains with significantly lower activity were Streptomyces spp., Bacillus megaterium, Bacillus aryabhattai, Acinetobacter spp. and Pseudomonas spp., while the Mitsuaria spp. and Rhizobium lusitanum strains displayed the lowest activity ( Figure 2 ). These last two strains were excluded from subsequent stages of the study. Rhizobium spp. emerged as members among the most active strains; however, this species is not a typical endophyte and was excluded from further stages of the study.

All experimental factors (3 phosphorus concentrations, 2 willow species and 2 bacterial inoculations) showed statistically significant effects on most plant growth parameters ( Table 1 ; Supplementary Figures 1, 2 ). A slightly better plant growth stimulation effect was observed for Paenibacillus spp. than for the other species (B2). Significant differences were also observed in response to different concentrations of phosphorus in the soil substrate. The positive or negative effects of the bacteria were dependent on the parameters tested, and the direction of change varied between the willow species and the bacterial variants (e.g., significant inoculation x species interactions listed in Table 1 ). High phosphorus concentrations in the soil significantly increased root fresh weight and length, while a decrease in leaf fresh weight was observed ( Table 1 ). Soil phosphorus effects on plant growth patterns varied between the two willow species (e.g., the significant species x factor interaction effects in Table 1 ). Changes in soil phosphorus concentration did not affect the phosphorus content in the leaves ( Supplementary Figure 3 ). Only B1 showed a significant decrease in leaf phosphorus content for the LP and HP treatments compared to the NP treatment for both Loden and Tora ( Supplementary Figures 3A, B ).

Regardless of experiment variant 55% of reads from Loden species and 70% for Tora was mapped to reference genome used in this study. Additionally, PCA analysis was performed to check for intergroup differences between samples ( Supplementary Figure 4 ). A transcriptome analysis of two willow species was performed to determine the effects of inoculation with the two selected PSB (B1 and B2) and different concentrations of phosphorus compounds (NP, LP and HP). The most distinctive difference was observed for plants inoculated with B2 of the Loden variety compared with the Ctr, where 2147 genes were downregulated and 2296 genes were upregulated. In addition, a similar trend was observed for both willow species in the HP and LP treatments: there was a greater number of genes that were upregulated compared to those in the NP treatment, the number of which ranged from 426 to 1920 genes ( Supplementary Figure 5 ). No significant trend was observed in terms of downregulation of genes. The average numbers of all genes in Loden and Tora whose transcription patterns differed from those of the control variant were approximately 2500 and 1000, respectively. According to Venn diagrams ( Figures 3A–F ) for the Loden species, a pool of genes common to both the Ctr and B1 and B2 variants was lower in the NP (15009 genes) treatment compared to the LP (18099 genes) and HP (18054 genes) treatments ( Figures 3A–C ). In addition, a very low number of common genes between the B1 variants and the Ctr were observed in the LP treatment for the Loden ( Figure 4B ) and Tora ( Figure 3E ) variants, which were 355 and 430, respectively, the lowest values recorded. As shown in Figures 4 – 6 , a Gene Ontology (GO) enrichment analysis was performed, and the 30 most differentially expressed genes associated with p values under 0.05 are presented for three phosphorus concentrations and two inoculation variants (bacteria B1 and B2 are compared to the corresponding non inoculated Ctr). Analyzing the effect of different phosphorus concentrations on the Loden willow species, we characterized the molecular functions in the NP treatments according to the presence of a large number of genes responsible for tetrapyrrole binding, haem binding, iron ion binding, GTP binding or ribonucleotide binding ( Figures 4A, B ). The LP treatment revealed a high number of genes responsible for the process of translation, and the HP treatment displayed a high number of genes associated with cytoskeleton formation and rearrangement ( Figures 5A, B , 6A, B ). In the case of the inoculation effect, the Loden variants inoculated with B1 (rather than B2) demonstrated a high expression of genes involved in cell membrane structure, phosphatase activity and acid phosphatase activity ( Figures 4 , 6 , 7 ). The variants inoculated with B2 were characterized by genes involved in hydrolase activity and pyrophosphatase activity and many genes related to cation and ion transport ( Figures 4 , 6 , 7 ). For the LP treatment, significant changes in the molecular functions of plants inoculated with B1 and B2 were observed, although some similarities were detected. In response to this variant, a significant increase in the number of genes responsible for DNA binding and regulation of transcription and translation was observed ( Figures 6A, B ). However, enrichment of GO terms related to direct phosphorus solubilization was not observed. Similar to the NP treatment, a number of genes responsible for the biosynthesis and metabolism of peptides that make up cell membranes were noted. In the category of biological processes, the variants inoculated with B1 under LP revealed processes related to the stress response (especially oxidative stress) and the plant defense system ( Figure 6A ). In the HP treatment, a shift from transcription to translation and, in particular, to ribosome synthesis was observed for the genes ( Figures 6A, B ). Moreover, genes related to chromatin formation and the formation of ribonucleoprotein complexes were found. For the NP, LP and HP treatments, GO terms associated with cell membranes were detected among the differentially expressed genes. The plants inoculated with B2 for HP were also characterized by genes related to different types of hydrolases, which were not observed in response to B1.

For the Tora willow species under NP, a high frequency of genes responsible for transcription, DNA binding and organic acid metabolism was observed ( Figures 4C, D ). Plants under the LP treatments were characterized by genes associated with transport and pollination, whereas under HP, a high frequency of genes associated with chloroplasts, thylakoids and photosystems I and II was observed ( Figures 6C, D , 7C, D ). For the effect of inoculation, B1-inoculated plants showed molecular functions belonging to tetrapyrrole and haem-binding activity and biological processes related to pollen recognition and formation and intercellular communication ( Figures 4C , 6C , 7C ). The plant inoculated with B2 was mainly characterized as expressing genes involved in metabolic processes and the biosynthesis of carboxylic acids, oxyacids and organic acids ( Figures 4D , 6D , 7D ). The Tora willow species under LP showed a change in gene expression and significant differences between the B1 and B2 variants ( Figures 6C, D ). In the case of the B1 variant, a high number of genes responsible for the stimulus response even at the cellular level was noted. A large number of genes responsible for transport also was detected. B2 under LP elicited the expression of a large number of genes responsible for transcription ( Figure 6D ). The he expression of genes related to pollen production, pollination and modification of various proteins was also evident, albeit with lower numbers. As with the LP treatment, little similarity was observed between B1 and B2 inoculation (with the exception of genes encoding transferases) ( Figures 6C, D ). Inoculation of plants with B1 resulted in high expression of genes responsible for transferase activity and coenzyme binding. A large number of genes related to the metabolic activity of various acids, glucan, lipids and carbohydrates also emerged. Inoculation with B2 increased the expression of genes involved in photosynthesis (thylakoid structure, photosystems I and II) and the synthesis of structural elements of plant cells. For each of the NP, LP and HP treatments, upregulation of a large number of genes responsible for the activity of different types of transferases was observed.

Among the 10 most frequently up- or downregulated genes found in each experimental variant, a general trend based on the willow species/species was observed in which genes for the Loden willow species were downregulated and genes for the Tora willow species were upregulated regardless of the experimental conditions tested (bacterial inoculant B1 or B2 and phosphorus concentration). This was particularly evident for genes associated with biological processes ( Figure 8 ) and molecular functions ( Figure 9 ).

Analysis of GO biological processes ( Figure 8 ) revealed the occurrence of 86 different up- or downregulated genes across all factors tested in the experiment. Among these genes, 43 were specific only to Loden, and 37 were specific to Tora. There were also 6 genes that were either up- or downregulated in both plant species tested. These genes were responsible for the cellular carbohydrate metabolic process, glucan metabolic process, glucan metabolic process, defense response, response to auxin, and lipid biosynthetic process. In response to both inoculation variants (B1 and B2), 13 common genes were observed, mainly those related to the metabolism and biosynthesis of fatty acids, glucan, and carboxylic acids. A strong effect of phosphorus (NP, LP, HP) was observed for fatty acid and monocarboxylic acid biosynthesis and metabolism, for which downregulation was observed in the NP treatment and upregulation in LP and HP treatments.

In terms of the GO cellular components, 65 different genes were observed, of which 28 were specific only to Loden (mainly downregulation), and 27 were specific to Tora (mainly upregulation) ( Figure 5 ). At 10, the pool of common genes for cellular components was larger than that for biological processes and included a large number of genes related to cellular elements involved in photosynthesis, i.e., thylakoids, photosystems I and II, the photosystem II oxygen-evolving complex and photosynthetic membranes. In response to different inoculants (B1 vs. B2), we observed the occurrence of 16 common genes, which were associated with terms related to chromosomes, the apoplast, DNA packaging or extracellular elements in addition to photosystem elements. For genes associated with cellular components, we also observed a strong effect of phosphorus on photosystems I and II, thylakoids, the photosystem II oxygen-evolving complex, and the oxidoreductase complex; moreover, the NP treatment downregulated the expression of these, in contrast to the LP and HP treatments, which upregulated the expression of these genes.

For GO molecular functions, our analysis revealed the presence of 73 genes, among which 37 were expressed only in Loden, and 27 were expression in Tora ( Figure 9 ). The pool of common genes comprised 11 genes and included those mainly responsible for calmodulin, haem and calcium binding, transferase activity, transcription and DNA-binding processes. We observed that the gene pool shared by Loden and Tora was characterized by the lack of an apparent division between downregulated and upregulated genes. The pool of common genes that responded to both B1 and B2 totaled 12. The abovementioned group of genes was joined by microtubule-binding activity, DNA-binding transcription factor activity and copper ion binding. In the case of molecular functions, the influence of phosphorus on gene regulation was observed only in the case of coenzyme binding, where the HP treatment elicited upregulation, which contrasted with the NP treatment, which caused those genes to be downregulated. Additionally, downregulation of two genes related to phosphatase activity was characteristic for the Loden B1 NP treatment combination.

From a pool of 64 PSB, two strains with the highest TCP-solubilization potential were selected (based on two screening stages carried out on solid and liquid TCP media) to investigate the effects of inoculation on the growth of willow seedlings under different phosphorus concentrations. The willow plants were grown in a substrate supplemented with TCP (Ca³(PO⁴)²), which is not a readily available source of P, and the phosphate solubilization process carried out by the bacterium inoculants was the only source of available P, with the exceptions of strains Bacillus aryabhattai (GL2_5_ED2) and Caballeronia glathei (GTL2_4_RP10), which solubilized TCP and DCP equally well. Bacillus spp. are well-known and well-documented PSM and have been tested many times for their use as potential biofertilizers of crop plants (Tahir et al., 2018; Ahmad et al., 2019; Miljaković et al., 2020; Mosela et al., 2022). In the present experiment, among the 15 most effective strains that solubilized phosphates, Bacillus spp. also were predominant (7 out of 15). This study found that Rhizobium sp. bacteria solubilized phosphate compounds most effectively, which highlights their great potential not only for nitrogen fixation but also as free-living soil bacteria that can solubilize P. Recent literature has confirmed the high efficiency of Rhizobium applied to crop plants, but there are no reports of the role of this species in tree growth (Afzal and Bano, 2008; Jaybhay et al., 2017; Verma et al., 2020; Mir et al., 2021; Shome et al., 2022).

In the pot experiment, strains Pantoea agglomerans (B1) and Paenibacillus sp. (B2) that effectively solubilized phosphates in both LP and HP media (analyzed by the molybdenum blue method) were selected. Pantoea agglomerans is described as a plant growth-promoting rhizobacterium (PGPR) capable of auxin biosynthesis, ACC deaminase production, ammonia production, phosphorus solubilization and increasing plant tolerance to salt stress (Majumdar and Chakraborty, 2015; Cherif-Silini et al., 2019; Emami et al., 2021). Because willow trees do not constitute a food crop, the potential pathogenicity of willow to humans was not an issue (Mackiewicz et al., 2016). The second strain we selected was Paenibacillus sp., which is described in the literature as a bacterium capable of nitrogen fixation; phosphorus solubilization; stimulation of the plant defense system; and production of ammonia, HCN, IAA and siderophores (Kumari & Thakur, 2018; Hussain et al., 2020). Both bacterial species have been previously experimented with as bioinoculants of plant species such as barley (Canbolat et al., 2006), wheat (Hussain et al., 2020), sugarcane (Quecine et al., 2012), poplar (Vaitiekūnaitė et al., 2021) and palm trees (Saadouli et al., 2021).

Inoculation with PSB and phosphorus treatments had significant effects on the growth patterns of the tested plants. However, it was difficult to determine whether the changes were caused by different soil phosphorus levels, bacterial inoculation, or a combination of both. Growth responses were different for Loden and Tora. For Loden, an increase in shoot thickness was found in response to higher soil P, while for Tora, the shoot length increased. In the NP treatment, the decrease in the leaf fresh weight could have been the result of phosphorus deficiency, because phosphorus deficiency has been shown to stimulate root allocation at the cost of leaf allocation (Ericsson, 1995). However, changes in phosphorus allocation to the roots or the leaves could also be related to changes in substrate pH, because increasing TCP concentrations may have increased the substrate pH. According to Walle et al. (2007), the optimum pH for willow growth is between 5 and 7.5, and the process of phosphorus solubilization by microorganisms is often associated with soil acidification (Liang et al., 2013). While alkaline phosphorus compounds are being formed, the pH becomes equilibrated by the process of microbial phosphate solubilization (Khan et al., 2014).

Plant transcriptome studies are becoming an integral part of all analyses aimed at determining the impact of biotic or abiotic stresses on various biological processes, cellular components or molecular functions. There are a growing number of publications examining the transcriptomic responses of tree crops and the impact of factors such as drought (Pucholt et al., 2015; Jia et al., 2020; Xu et al., 2021), salinity (Yao et al., 2018; Zhou et al., 2020; Pang et al., 2022) and pest resistance (Wang et al., 2020). To date, however, there are no publications showing differences in gene expression under P-deficient conditions. Our results allowed us to determine the effects of phosphate-solubilizing bacterium inoculation with increasing concentrations of plant-available phosphate on the expression of genes related to the phosphorus stress response and phosphorus metabolism, mobilization and remobilization.

Our transcriptomic analysis showed high variability between the gene expression in the leaves of the two willow species tested (Loden vs. Tora), the inoculation variants (B1 and B2) and the phosphorus concentrations in the substrate (NP, LP, HP). The results showed that Loden exhibited approximately 1500 more differentially expressed genes than did the Tora willow species. As our previous study has shown, these species differ significantly in their morphology (Hoeber et al., 2017) and root-associated microbiome (Koczorski et al., 2021a; Koczorski et al., 2022). As reported by Hoeber et al. (2017), Tora is characterized by highly efficient biomass production (tall shoots) but a smaller leaf area. Conversely, Loden is characterized by a large leaf area but relatively low amount of biomass production. Moreover, under NP conditions, Loden was shown to express a low number of genes in common between both the B1 and B2 variants, suggesting that inoculation significantly affects gene activity, specifically under these conditions. According to previous studies, one of the plant responses to inoculation is increased activity of genes responsible for DNA replication and cell division (Camilios-Neto et al., 2014). In addition, the activity of genes related to response to stimuli, signaling, regulation of biological processes and metabolic processes were observed (Liu et al., 2012). In our analysis, similar gene activities in both the Loden and the Tora willow species, especially under LP conditions, were observed. Additionally, the expression of genes responsible for photosynthesis, transcription and translation in the LP and HP treatments was also observed. This was manifested by the activity of genes such as those involved in the glucan metabolic process, thylakoids, photosystems I and II, translation, DNA-binding translation factor activity, transcription regulation and sequence-specific DNA binding.

Phosphorus plays a very important role in the photosynthetic carbon reduction cycle, also known as the Calvin-Benson cycle, because the final product of this process is 3-phosphoglyceraldehyde (triose-P), which is then (already in the form of sucrose) transported to the cytosol of the cells or remains in the chloroplasts where it is converted into starch (Rychter and Rao, 2005). Among the previously mentioned genes are also those responsible for glucan metabolic processes. Glucan phosphatase plays an important role in the transient metabolism of starch in leaves, as glucan phosphatase is the main storage material and consists mainly of glucose polymers. During the diurnal cycle of photosynthesis, starch is deposited during the day to be used up at night to provide sufficient metabolites needed for plant growth (Meekins et al., 2016). Our results suggest that the bacteria present in the substrate significantly improved the availability of starch and thus had a positive effect on plant growth. For the Tora willow species inoculated with B1, activity of genes responsible for pollen recognition and formation (NP and LP treatments) was observed. The direct effect of inoculation on pollen formation has not yet been investigated, but according to a study by Lankinen et al. (2000), phosphorus has a positive effect on the quantity and quality of pollen formed, while abiotic stresses such as salinity do not affect this process. Among the variants tested, only the Loden willow species inoculated with B1 under LP showed the expression of genes directly related to stress and defense mechanisms. Studies conducted on Arabidopsis thaliana leaves under P-deficiency conditions showed that these plants exhibit high expression of genes involved in phosphatase synthesis and genes that encode ribonucleases and sulfolipid biosynthesis proteins that replace P, contributing to phosphorus remobilization (Müller et al., 2007). Similarly, we detected genes involved in ribonuclease activity and phosphoric ester hydrolase activity, which may indicate P-deficiency stress. The genes indicated above (except for those encoding ribonucleases) were also expressed at lower numbers in response to the other treatment combinations in this study, which may be related to the rapid growth of the plants at an early stage of development. As a woody plant species, willows are known for their rapid growth, and as research has indicated, phosphorus remobilization and phospholipid replacement by sulfur and galactolipids is a strategy employed to increase photosynthetic efficiency (Lambers et al., 2012).

The information obtained in this analysis largely overlaps with the results obtained from the GO enrichment analysis, providing additional information about genes that are upregulated and downregulated. The first and most prominent result of the analysis is the observed trend, in which the genes of the Loden willow species were downregulated and those of the Tora willow species were upregulated. As mentioned earlier, this may be a result of the difference in biomass allocation between the two willow species (Ericsson, 1995).

For the Tora willow species, a large effect of phosphorus on genes related to photosynthesis, thylakoids and starch metabolism in the form of upregulation of photosynthesis, photosystems I and II, thylakoid parts, carboxylic acid, carbohydrates, glucan metabolism and biosynthesis genes was found. All of the abovementioned genes are highly dependent on phosphorus because it is a component of ATP, NADPH, nucleic acids, sugar phosphates and phospholipids (Carstensen et al., 2018). The effect of phosphorus on carboxylic acids seems to be substantial due to the important role of these compounds as plant signalling molecules; also important are compounds like malate, which can be converted into NADH and NADPH, thus increasing the efficiency of photosynthesis (Mutz et al., 2013). Additionally, the presence of carboxylic acids has been confirmed in the apoplast, where activity related to the regulation of stomatal opening occurs (Meyer et al., 2010).

Among the downregulated genes for the Loden willow species, genes involved in lipid peptide metabolism and biosynthesis, nucleic acids, and ion and metabolite transfer were detected, while upregulation of genes involved in the defense mechanism and response to stress was also observed. Gene expression studies conducted on plants subjected to P-deficiency stress have indicated that the most common plant response to phosphorus deficiency is the downregulation of genes involved in ATP synthesis, translation, transport, carbohydrate synthesis and photosynthesis, which was observed for all the treatment combinations involving Loden (Bai et al., 2014; Chu et al., 2018). In addition, the upregulation of genes related to sulfur compounds was observed for the B2-inoculated Loden treatment combinations, indicating an attempt to replace phosphorus to increase phosphorus-use efficiency (Lambers et al., 2012). Notably, plants in the LP and HP treatments did not show signs of phosphorus deficiency, so the phenomena are related to species differences and biomass allocation rather than phosphorus itself.