Deciphering the Structural and Functional Diversity of Rhizobacteria from Stone Pine Inoculated with Plant Growth Promoting Rhizobacteria (PGPR) before and after Transplanted into Degraded Agricultural Soil

: The use of plant growth-promoting rhizobacteria (PGPR) inoculated on plants has shown that it can increase the success of reforestation and accelerate soil recovery by improving soil microbial diversity. Three PGPR isolated from natural pine populations were selected for their metabolic capabilities and taxonomic affiliation (Z4.3; Bacillus sp., Z5.4; Arthobacter sp., and Z7.15; and Pseudomonas sp.) when inoculated alone or in combination (consortium) on stone pine seedlings before transplanting to the field. Before transplanting and after nine months, rhizospheric soil samples were collected for structural and functional metagenomic studies. First, the data were analyzed using EasyMAP. Neither alpha nor beta diversity showed significant differences between the samples, although unique taxa representative of each sample were detected. The predominant phylum in all cases was Proteobacteria, followed by Bacteroidetes and Acidobacteria. The linear discriminant analysis (LDA) effect size (LEfSe) found significantly over-represented taxa in some samples, highlighting different representatives of the order Sphingomonadales in several of them. Functional inference performed with PICRUSt also showed significantly over-represented functions in some samples. The study demonstrates that PGPR have a positive effect on plants and cause detectable changes in microbial communities in terms of both structure and function.


Introduction
Soil degradation is probably one of the most important environmental problems we face, and it is estimated to have negative consequences for the well-being of at least 3.2 billion people.Agricultural activities, particularly the abandonment of agriculture, are important drivers of soil erosion and degradation [1][2][3].Climate change is very likely to intensify the process of soil abandonment; so, proper environmental management is essential to prevent further soil degradation.A very important effect of soil degradation is the loss of biodiversity.Erosion has a negative effect on microbial diversity and soil functionality [2].Soil microorganisms are essential for the functioning of ecosystems since they regulate energy flow and mass flows and drive soil ecosystem responses to anthropogenic disturbances and environmental changes.Any change in microbial diversity has a significant environmental impact [4][5][6][7].
The reforestation of degraded soils is being widely used as a restoration strategy for this type of soil; it is a strategy that can have a very positive influence on the improvement of the functioning of the ecosystem and the physical-chemical properties of the soil, but it can also decisively influence the dynamics of the soil microbial community [8].This approach is not free of difficulties; one of the most important difficulties is related to the problems the transplanted forest species face in trying to colonize these soils, which are sometimes severely degraded and have lost some of the qualities required for the plants to root.
Root-colonizing bacteria, known as rhizobacteria, have developed symbiotic relationships with nearly every type of plant.Often known as plant growth-promoting rhizobacteria (PGPR) [9], these bacteria can increase plant growth and productivity with a variety of mechanisms, including the synthesis of siderophores and phytohormones, the solubility of inorganic phosphate, and antifungal activity.In general, when plants grow in poor and/or stressed soils, the positive impacts of PGPR on plant development and health are more noticeable [10].The ability of PGPR to colonize, survive, and remain in the rhizosphere, together with their interactions with the local microbial community, are widely acknowledged to be key factors in their success.The impact of PGPR inoculation on the native microbial community and activity in either bulk or rhizospheric soil has therefore received increasing attention [11].
Planting forest tree seedlings treated with beneficial microorganisms that promote plant growth and soil health is one of the techniques advised for the reclamation and restoration of wasteland [12].Beneficial microbe-inoculated trees are typically chosen for planting on damaged land because of their capacity to withstand biotic and abiotic challenges and boost soil fertility [13][14][15].By increasing the nutrient status and quality of the soil, these advantageous microorganisms that form colonies in and around root systems support plant growth.Research has demonstrated the favorable effects of applying microbial consortia made up of particular microorganisms on plant growth, particularly when the microorganisms are inoculated in nursery beds or root trainers during the seedling stage [16,17].
The stone pine (Pinus pinea L.), a distinctive species of Mediterranean flora, has been the subject of study.Its estimated global size is 600,000 hectares, of which more than 400,000 Ha are in Spain.Because of this species' economic significance, primarily in the production of pine nuts and timber, it has been employed since antiquity [18].Because its roots are usually well developed, the super-producing stone pine plays a crucial role both in fixing coastal dunes and in controlling erosion in mountainous regions.This pine species is perfect for recovery efforts for degraded soils because of all these qualities [19].
Based on the above, a screening of bacteria associated with the stone pine rhizosphere of populations already established in degraded soils for 5 years was carried out.As a result of this screening, 268 isolates were obtained, and morphological (Gram-positive sporulated, Gram-positive non-sporulated bacilli, and Gram-negative cocci) and in vitro PGPR traits were evaluated.Forty strains were selected to sequence the 16S rRNA gene and to determine the strains with the highest PGPR capacities from all the zones and parataxonomic groups.Three strains were chosen to inoculate one-year-old pine trees that were going to be transplanted to abandoned agricultural land.
The main objective was to determine a procedure to improve plant adaptation to degraded agricultural soils and to improve soil biodiversity by delivering PGPR inocula to seedlings during the nursery process.To achieve this purpose, a screening of PGPR from stone pine growing in low-nutrient soils was carried out; the isolates were characterized according to genetic diversity and putative beneficial traits.Selected strains were delivered to stone pine seedlings in the nursery.After transplanting to agricultural soil, the effects of the inoculum on the rhizosphere from a structural and functional perspective and on plant growth were determined.

Experimental Design
Bacteria were isolated from the rhizosphere of Pinus pinea from ten different zones in Segovia, Spain (41 • 14 ′ 29.2 ′′ N 3 • 50 ′ 18.2 ′′ W).All the isolates were characterized in parataxonomic groups (Gram-positive non-sporulated bacilli, Gram-positive sporulated bacilli, Gram-negative bacilli, and Gram-positive cocci), and in vitro putative PGPR traits tests were performed.From each zone and each parataxonomic group, the strains with the highest PGPR capacities (qualitatively and quantitatively) were chosen to sequence the 16S rRNA gene.Based on PGPR capabilities, while trying to eliminate genetic redundancy, three bacteria of three different genera were chosen: Z4.3, Z5.4, and Z7.15.
The plants were inoculated six times, from April to September 2019, once a month.During this period in the nursery, the height and diameter of the base of the stem were measured in April (before starting inoculations) and in June.
The plants were transplanted to the field in October 2019.Rhizosphere samples were taken in the nursery before transplanting in September 2019 (autumn samples) and after transplanting the trees to the degraded soil in June 2020 (spring samples); the samples were taken from the rhizosphere of the transplanted pines to carry out a diversity analysis, and the height and diameter were measured again.

Isolation of Bacterial Strains from the Rhizosphere of Pinus pinea
In October 2018, pine rhizosphere soil was collected from ten zones on land belonging to "El Ejidillo Viveros Integrales" (41 • 14 ′ 14.6 ′′ N 3 • 50 ′ 49.3 ′′ W).Zones 1 to 6 were in a plot (8 Ha) that was reforested with stone pine 14 years ago.Zones 7 and 8 had had pine trees for 8 years.Zones 9 and 19 had had pine trees for 6 years.
The rhizospheric soil (soil deeply attached roots) of four plants was randomly pooled in plastic bags and transported at 4 • C to the laboratory.In an omnimixer, 2 g of rhizosphere soil were suspended in 2 mL of sterile distilled water and homogenized for 1 min.Five hundred µL of the soil suspension was plated on a plate count agar, PCA (Condalab, Madrid, Spain), and incubated for 48 h at 28 • C. One hundred µL of the soil suspension was used to prepare successive ten-fold dilutions in a final volume of one mL.From each zone, 50 distinct colonies were isolated.Each colony was separated and categorized into four parataxonomic categories based on physical traits, sporulating capacity, and Gram staining: Gram-positive endospore-forming bacilli, Gram-positive non-endospore-forming bacilli, Gram-negative bacilli, and Gram-positive cocci.All the isolates were kept at −20 • C in nutrient broth with 20% glycerol.

PGPR Potential Traits Tests
All the isolates were tested for in vitro production of siderophores, chitinases, auxins, and phosphate solubilization.
Using the medium outlined by Alexander and Zuberer [20], siderophore synthesis was detected.In this medium, iron is bound to CAS (Chrome azurol S) as a Fe-CAS complex.The iron from the Fe-CAS complex is separated by bacteria that can produce siderophores.After 168 h of incubation, if the isolate was able to produce at least 5 mm of a yellow to orange colored halo surrounding the colony, it was deemed positive for siderophore production.
The strains were grown on a culture medium containing colloidal chitin, K 2 HPO 4 , KH 2 PO 4 , MgSO 4 , NaCl, KCl, yeast extract, and agar in order to evaluate the strains' capacity to generate chitinases [21,22].After 5 days of incubation at 30 • C, the bacteria able to synthetize chitinases produced a transparent halo (the medium is opaque) around the bacterial growth zone.This halo was measured in millimeters.After 5 days of incubation, if the isolate could form a transparent halo around the colony of at least 5 mm, it was deemed positive for chitinase production.
A colorimetric approach was employed for the identification of auxin-like compounds [23].The bacterial isolates were plated in a grid arrangement on half-strength Tryptic Soy Agar (TSA; Difco Laboratories, Sparks, MD, USA) and were allowed to grow for Soil Syst.2024, 8, 39 4 of 19 24 h.Then, a nitrocellulose membrane with a diameter of 82 mm (Amersham Biosciences, Little Chalfont, UK) was placed over the inoculated bacteria in the plate.Salkowski's reagent (2% (v/v) 0.5 M FeCl 3 in 35% (v/v) perchloric acid) was then applied to the membrane, and it was left to develop color for 2 h.Using known quantities of pure IAA, the assay's sensitivity was ascertained (Sigma Chemical Co., St. Louis, MO, USA).
The ability of the bacteria to solubilize phosphate was determined by plating the bacteria in potato dextrose yeast extract agar (PDYA, pH 7.0) medium with freshly precipitated CaHPO 4 [24].To produce a precipitate of CaHPO 4 , 50 mL of K 2 HPO 4 10% (w/v) and 100 mL of CaCl 2 10% (w/v) were added per liter of sterile PDYA.After a 24 h incubation period at 28 • C, the hydrolysis halo formed in this culture media around the colonies was measured in millimeters.The capacity of the strain to solubilize phosphate was confirmed by the existence of a hydrolysis halo.If the isolate could create a halo around the colony that measured at least two millimeters, it was deemed positive for phosphate solubilization.

16S rRNA Gene Sequencing
From each zone and each parataxonomic group, the strains with the highest PGPR capacities (qualitatively and quantitatively) were chosen to sequence the 16S rRNA gene; 40 strains were chosen in total.
Using the Ultraclean Microbial DNA Isolation Kit (Mo Bio, Carlsbad, California, USA), the DNA of every isolated strain was extracted.16S rRNA universal primers, 1492R (5 ′ TACGGYTACCTTGTTACGACTT3 ′ ) and 27F (5 ′ AGAGTTTGATCMTGGCTCAG 3 ′ ), were used to amplify each DNA sample.In the amplification reactions, 0.5 µL of primer F (30 µM) and 0.5 µL of primer R (30 µM), 2.5 µL of 10X standard reaction buffer with MgCl2 (Biotools, Madrid, Spain), 0.625 µL of dNTPs (10 mM each, Biotools, Madrid, Spain), 0.375 µL of 100% DMSO (dimethyl sulfoxide), and ultrapure water up to a volume of 25 µL were used.The reactions were first incubated for 2 min at 94 • C in a thermocycler (Applied Biosystems Gene Amp PCR system 2700).Following this, the reactions were exposed to 30 cycles, which included 94 • C for 30 s, 55 • C for annealing, and 72 • C for 1 min (with each cycle rising by 5 s from cycle 11 to cycle 30).After a final 7 min of incubation at 72 • C, the mixtures were sequenced in an ABI PRIMS ® 377 DNA sequencer (Applied Biosystems, Foster City, CA, USA), and the PCR products were purified using a QIAquick PCR purification kit from Qiagen.Clone Manager Professional Suite v6.0 was used for editing, while Sequence Scanner v1.0 was used to see the sequences.On the server MAFFT v6.0 (http://mafft.cbrc.jp/alignment/software/. Accessed on 15 January 2019), sequence alignment was performed and assessed by BLASTN 2.2.6 in the National Center for Biotechnology Information (NCBI: http://www.ncbi.nlm.nih.gov/.Accessed on 15 January 2019) and the Ribosomal Database Project Release 10 (RDP: http://rdp.cme.msu.edu/.Accessed on 15 January 2019) databases.

Phylogenetic Tree
An unrooted phylogenetic tree was constructed with MEGA v4.0.2, with the sequences aligned in MAFFT v6.0.Neighbor joining was used to infer the evolutionary distances.The evolutionary history of the taxa under study was assumed to be represented by the bootstrap consensus tree that was generated from 1000 replicates.Next to the branches are the proportion of replicate trees where the related taxa clustered together in more than half of the 1000 bootstrap replicates.Every position in the dataset with gaps and missing data was removed (full deletion option).The utilized sequences had accession numbers from PP033436 to PP033475 and were deposited in GenBank (NCBI).

Selection of Bacteria and Verification of the Viability of the Consortium
Based on PGPR capabilities and the attempt to eliminate genetic redundancy, three bacteria were chosen: Z7.15 (Pseudomonas sp.), Z4.3 (Bacillus megaterium), and Z5.4 (Arthrobacter sp.).
To check the compatibility of the three selected bacteria and to rule out the possibility that any of them inhibited the growth of any of the others, 3 equidistant drops of 10 µL from the 3 strains grown in nutrient broth for 24 h at 28 • C were placed in the center of a plate count agar (PCA, Condalab, Spain) plate.The plate was incubated up to 48 h at 28 • C to allow bacterial growth or inhibition halos to appear.

Inoculation of Selected Bacteria and Consortium in Nursery Pine Plants and Transplanted Plants
The plants were inoculated at root level with 100 mL of each bacterial culture diluted in water to 10 7 cfu/mL.For the consortium, the three bacteria were grown independently and mixed in equal proportions to reach 10 7 cfu/mL.The control plants were mock inoculated with 100 mL of water.
The inoculations were applied once a month between April and September 2019, making a total of 6 inoculations with 100 mL of 10 7 cfu/mL in each one.
During this period in the nursery, the height and diameter of the root neck were measured in April 2019 (before starting inoculations) and in June 2019.
In October 2019, the plants were moved to the field, where they were transplanted following a randomized block design.The field was very close to the nursery, where the plants had been in alveolus trays.
The plot where the pine seedlings were transplanted is owned by El Ejidillo Viveros Integrales (41 • 14 ′ 40 ′′ N, 3 • 50 ′ 10 ′′ W; 931 m).It is an abandoned agricultural plot that is not cultivated but is ploughed annually to keep it fallow.In the period between the two ploughing processes, a process of secondary ecological succession begins.The species that colonize the plot during this period include Papaver rhoeas, Portulaca oleracea, Erodium sp., Valerianella sp., Xanthium spinosum, Helianthus annuus, Hordeum vulgare, etc.According to the USDA classification, the soil has a loamy sand texture (sand: 74.31% ± 1.46; silt: 13.55% ± 1.18; clay: 12.14 ± 0.51).It has a pH of 7.99 ± 0.03, a nitrogen percentage of 0.125% ± 0.006, an organic matter percentage of 1.58% ± 0.07, and a C/N ratio of 7.72 ± 0.15.The climatology of the area, according to the data collected by the Cuellar weather station located a few kilometers away, has an average temperature of 12-13 • C and an annual rainfall of 450-500 mm.The climate of the area is continentalized Mediterranean, with very hot summers and very cold winters with temperatures below zero.

Analysis of Biological and Functional Diversity of the Rhizosphere
A metagenomic analysis of the 16S rRNA gene of the rhizospheres of all the treatments was performed with the samples before transplanting (nursery, autumn sampling time) and after transplanting (spring).In both sampling times, the rhizosphere samples of three plants of each replicate were taken.The samples were transported to the laboratory at 4 • C.
Firstly, DNA was extracted with the DNeasy PowerSoil DNA isolation kit.The DNA was sent at 4 • C to the Allgenetics company (La Coruña, Spain), which performed the metabarcoding libraries for the bacteria, sequencing using Illumina (Run 1/4 Illumina MiSeq PE300), and bioinformatic analysis to obtain the biodiversity of the bacteria.
EasyMAP, a user-friendly online platform for evaluating 16S ribosomal DNA sequencing data, was utilized to evaluate the 16S rRNA sequences [25].EasyMAP integrates modules from Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt), linear discriminant analysis effect size (LefSe) [26], and QIIME2.Quality control is the initial stage in EasyMAP, and the QIIME2 plug-in q2-dada2 was used for this purpose.In this stage, the maximum expected error (-p-max-ee) is two, and the maximum base number for each sequence is less than 1000.(-p-trunc-q) is the truncated quality value, and it is two."Consensus" is the chimera filtering method (-p-chimeramethod).A sequence being checked for chimericity (-p-min-fold-parent-over-abundance) must have a minimum of one possible parent.One million reads were used to learn the error model (-p-n-reads-learn).Next, the creation of phylogenetic trees and the taxonomic analyses were carried out."q2-feature-classifier" is the taxonomic analysis module that is employed, and "q2-alignment" and "q2phylogeny" are the phylogenetic tree-generating modules.In this stage, EasyMAP gives users the choice to eliminate non-target genes from the OTU table and representative sequences, such as sequences for mitochondria and chloroplasts.For alpha diversity analysis, EasyMAP produces "observed OTUs", "Faith's Phylogenetic Diversity", "Shannon index", and "Pielou's evenness".Microbes with significant differences across groups are identified using LEfSe, which EasyMAP applies to provide taxonomic differential abundance analysis.When conducting a Wilcoxon test between distinct subgroups, the default alpha value is 0.05.Two is the standard threshold for the logarithmic linear discriminant analysis (LDA) score's absolute value.For the Kruskal-Wallis test between statistical groups, the default alpha value is 0.05.Finally, the Greengenes database (version 13_8, 99% OTU) and VSEARCH [27] are used for close-reference clustering to predict the function of the microbes.This ensures that the representative sequences map to the Greengenes database.PICRUSt locates the matching function in the KEGG database using the mapped Greengenes ID.

Statistical Analysis
To examine the statistical differences in the stem diameter and height between the treatments in each sample period (June 19 and June 20), one-way ANOVA with repetitions was carried out.Statgraphics Plus 5.1 for Windows was used to verify the homoscedasticity and normality of the variance prior to ANOVA analysis, ensuring that all the analytical conditions were met.A Fisher test was employed when significant differences (p < 0.05) were observed.The sequences obtained in the metagenomic analysis were analyzed using an online platform for analyzing 16S ribosomal DNA sequencing data named EasyMAP (http://easymap.cgm.ntu.edu.tw.Accessed on 10 December 2020).

Results
Although the plan was to collect 50 bacteria per zone, this number was not reached because fewer bacteria grew.All those that grew on the plates were collected, for a total of 268 strains.The minimum numbers of isolates were collected from zones 6 and 4, where 6 and 11 strains were collected, respectively, and the maximum number of isolates was obtained from zone 9, with 48 isolates.The most representative parataxonomic group was the Gram-negative bacilli, followed by non-sporulating Gram-positive bacilli.
The predominant activity in the pine rhizosphere was the production of siderophores, followed by the solubilization of phosphates (Figure 1), both of which are related to the mobilization of nutrients.
Soil Syst.2024, 8, x FOR PEER REVIEW 7 of 18 annotation (S: siderophores production, P: phosphate solubilization, A: auxin production, and Q: chitinases production).Strain Z4.3 was a Bacillus megaterium that produced siderophores, solubilized phosphates, and produced auxins in much higher concentration than the other B. megaterium strain with these three capacities.Strain Z5.4 was the only strain of Arthrobacter sp.which, in addition to producing siderophores, produced auxins.Lastly, strain Z7.15 was the only Pseudomonas strain that combined three capacities (siderophore production, phosphate solubilization, and auxin production).At the beginning of the study, in April 2019, before starting the inoculations, the trees had an average height of 54.760 cm ± 0.903 and an average diameter of 1.107 cm ± 0.017, which was measured at the base of the stem.These parameters were measured again after 2 months (and 2 inoculations) in June 2019, and in June 2020 (6 inoculations and 9 months from transplanting), these parameters were measured again (Table 1).It is striking that we only detected auxin production in zones 8, 9, and 10, which corresponded to the nursery area, and the production of chitinases was only in four areas (Figure 1).
Based on the best performance in the PGPR tests in each sampling area, forty one strains were selected to sequence the 16S rRNA gene.The results were annotated in the NCBI database, and a phylogenetic tree was made with the sequences (Figure 2).Three main groups were found at the genus level: Pseudomonas, Bacillus, and Arthrobacter.One isolate was selected from each of these groups based on their PGPR capabilities.In the phylogenetic tree, the PGPR capacities of each strain are represented next to the annotation (S: siderophores production, P: phosphate solubilization, A: auxin production, and Q: chitinases production).Strain Z4.3 was a Bacillus megaterium that produced siderophores, solubilized phosphates, and produced auxins in much higher concentration than the other B. megaterium strain with these three capacities.Strain Z5.4 was the only strain of Arthrobacter sp.which, in addition to producing siderophores, produced auxins.Lastly, strain Z7.15 was the only Pseudomonas strain that combined three capacities (siderophore production, phosphate solubilization, and auxin production).
At the beginning of the study, in April 2019, before starting the inoculations, the trees had an average height of 54.760 cm ± 0.903 and an average diameter of 1.107 cm ± 0.017, which was measured at the base of the stem.These parameters were measured again after 2 months (and 2 inoculations) in June 2019, and in June 2020 (6 inoculations and 9 months from transplanting), these parameters were measured again (Table 1).Inoculation with the consortium significantly increased the trees' heights in June 2019, but the stem diameters were not modified by the treatments (Table 1).One year later, in June 2020, trees' heights were significantly increased by inoculation with Z5.4,Z7.15 and the consortium, while the stem diameters were increased by Z7.15 and the consortium inoculation (Table 1).
The rhizobacterial biological diversity calculated using metagenomic analysis was not significantly modified by the different treatments (Figure 3a).However, there was a trend of lower Shannon index values in the spring (after transplanting) than in the autumn (nursery) samples (Figure 3b), except in Z4.3, which was maintained (Figure 3a).The rhizobacterial biological diversity calculated using metagenomic analysis was not significantly modified by the different treatments (Figure 3a).However, there was a trend of lower Shannon index values in the spring (after transplanting) than in the autumn (nursery) samples (Figure 3b), except in Z4.3, which was maintained (Figure 3a). Figure 2. Phylogenetic tree constructed with the 16S rRNA sequences.The evolutionary distances were inferred using the neighbor joining method.The bootstrap consensus tree inferred from 1000 replicates was taken to represent the evolutionary history of the taxa analyzed.Annotations of bacteria included in the tree were obtained from the NCBI (http://www.ncbi.nlm.nih.gov/.Accessed on 15 January 2019).PGPR capacities of each strain are represented next to the annotation (S: siderophores production, P: phosphate solubilization, A: auxin production, and Q: chitinases production).The red arrows indicate the PGPRs selected to inoculate the pine seedlings.---------------------AUT UMN SPRING The bacterial species that were unique to the different treatments in both sampling times (before and after transplanting) are presented in Figures S1 and S2, respectively.Before transplanting (in the nursery), the highest number of unique species was found in the control (98 species were present exclusively in this treatment), while in the sampling time after transplanting, this was observed in the consortium (107 unique species), followed by the control (85 unique species).Comparing each treatment in both sampling times, the number of unique species found was maintained, except in the Z4.3 and consortium treatments, in which the number of unique species was higher in the field (spring) than in the pot in the nursery (autumn).
The most represented phyla were Proteobacteria, Bacteroidetes, Acidobacteria, and Actinobacteria (Figure S3).In the Proteobacteria and Bacteroidetes phyla, the abundance was higher in the samples after transplanting than before and was outstanding in the control and Z7.15.It is also worth highlighting the abundance of taxa from the Firmicutes phylum in treatment Z4.3 in the spring sampling time.
Figures S4-S6 show the most abundant OTUs at the class and order level of the four phyla indicated above.It can be seen that at these taxonomic levels there is a clear relationship between the abundance of the different taxa and the treatment applied (different PGPR or consortium), as there are classes or orders that are more or less abundant in the microbiome of the plants studied, depending on the treatment.
Figure 4 shows the linear discriminant analysis (LDA) and shows which taxa are significantly more represented in the different treatments.In general, the different taxonomic groups corresponding to the Gram-negative bacteria were the most represented, except for the rhizosphere of the plants inoculated with Z4.3 before transplanting (autumn), where Gram-positive bacteria, especially Bacillus, were the most represented group.Among the Gram-negative bacteria, the taxa corresponding to the alpha-proteobacteriaceae, such as the Caulobacteraceae family or the Sphingomonadales, were abundant; the latter were especially abundant in several of the treatments (Z5.4 after transplanting, Z4.3 before transplanting, or the consortium after transplanting).Some of the significantly represented taxa were related to serious opportunistic infections of humans, such as the genera Legionella (Z7.15 in spring), Massilia (consortium before transplanting), and some species of the Sphingomonas genus (consortium after transplanting).In the rhizosphere of plants inoculated with Z7.15 after transplanting, the genus Chthoniobacter was significantly represented.
plants inoculated with Z7.15 after transplanting, the genus Chthoniobacter was significantly represented.Figure 5 shows the statistically more represented microbial functions in the different treatments inferred by PICRUSt from the 16S rRNA sequences.The treatment with the higher number of represented functions is Z7.15 before transplanting, including carbohydrate metabolism and secondary metabolite synthesis, among others.Functions related to Figure 5 shows the statistically more represented microbial functions in the different treatments inferred by PICRUSt from the 16S rRNA sequences.The treatment with the higher number of represented functions is Z7.15 before transplanting, including car-bohydrate metabolism and secondary metabolite synthesis, among others.Functions related to neurodegenerative diseases, such as Parkinson's disease (Z5.4 after transplanting), Alzheimer's disease (Z7.15 after transplanting), or Huntington's disease (control after transplanting) also appeared in several treatments.In the rhizosphere of plants inoculated with Z4.3 before transplanting, the functions related to important physiological processes in Bacillus sp., such as sporulation and germination, were represented.Cluster analysis of bacterial OTUs based on Jaccard distance (i.e., PCoA; principal component analysis) was used to evaluate the similarity of the bacterial community composition (beta diversity).First, a PCoA with the samples from the two sampling times was performed.A clear distribution of samples by sampling time was observed, but there was no grouping pattern within each group (data shown).For this reason, an independent PCoA was made for each sampling moment.The cumulative variance contribution of the Cluster analysis of bacterial OTUs based on Jaccard distance (i.e., PCoA; principal component analysis) was used to evaluate the similarity of the bacterial community composition (beta diversity).First, a PCoA with the samples from the two sampling times was performed.A clear distribution of samples by sampling time was observed, but there was no grouping pattern within each group (data shown).For this reason, an independent PCoA was made for each sampling moment.The cumulative variance contribution of the first two principal components extracted in the first sampling time (before transplanting, in autumn) was 16.83%, which accounted for 9.065% and 7.77% of the variance of the variables, respectively.In the second sampling time (after transplanting, in spring), the cumulative variance contribution of the two principal components was 25.54%, which explained 14.10% and 11.44% of the variance of the variables, respectively.The PERMANOVA performed did not indicate significant differences between any of the treatments; this was probably due to the low variation explained by the analyses: a group formed by the three control replicates, another group formed by the Z5.4 treatment replicates, and the last group formed by the rest of the samples (Z4.3, Z7.15, and the consortium).In the PCoA performed with the samples of the second sampling time (after transplanting, in spring), there were not such clear groupings, except for Z4.3 and Z7.15, which remained in the same cluster.

Discussion
Pine plants grafted with super-productive pine branches are one of the strategies used to increase pine nut production.However, transplanting these trees to degraded soil can present survival and production problems; so, growing plants with good fitness and good adaptation capacity is of enormous interest.The application of beneficial rhizobacteria is presented as a biotechnological tool to improve the adaptation and growth capacity of these plants with a relevant economic impact.On the other hand, reforestation of degraded agricultural land with these pine trees inoculated with PGPR can contribute to the recovery of soil biodiversity.At the same time, it is interesting to analyze the effect of inoculation on the diversity of soil microorganisms; this diversity is key to the health and fertility of the soil and plant growth [6,7].
In the phylogenetic tree constructed with the 16S gene sequences of the rhizobacteria (Figure 2), the strains isolated from different areas grouped together in the tree and showed 100% similarity; these included strain 33 (Z7.15), 39 (Z8.14), and 47 (Z9.43) in the Pseudomonas group or 28 (Z6.1) and 31 (Z7.4) in the Bacillus group.This indicates the capacity of plants to recruit strains that provide the best advantage for survival.This is not uncommon and has been described in the literature on multiple occasions [28][29][30][31].
A strikingly high percentage of the isolated pine rhizobacteria showed activities related to nutrient mobilization, like production of siderophores and solubilization of phosphates (Figure 1).The plant itself has a complex system for mitigating stresses brought on by degraded land; this system is further enhanced by many root-associated bacteria, which are able to trigger a variety of effective coping strategies simultaneously [32].The ability of plants to select rhizobacteria is determined by the type of exudates they exude from their roots.In addition to sugars, plants exude other selective molecules of different chemical natures to attract those microorganisms that are useful in certain circumstances [33][34][35].For example, Arabidopsis plants have been shown to exude coumarins to attract microorganisms capable of producing siderophores in situations of iron scarcity [36].The case of iron is not the only one; the selection of microorganisms that help the plant to mitigate nutrient stress situations can be performed for other nutrients [37,38].Many rhizosphere bacteria can enhance plant growth by supplying essential compounds like phosphorus [39], and this is achieved by several mechanisms [40], such the secretion of organic acid that lowers pH, the chelation of ions bound to P, and competition with P for adsorption sites in the soil or the production of enzymes to solubilize P [41].
After two months and two inoculations, the consortium significantly increased the trees' heights (Table 1).After one year (after six inoculations), the heights were significantly increased by Z5.4,Z7.15, and the consortium, and the trees' stem diameters were increased by Z7.15 and the consortium (Table 1).Therefore, the selected PGPR strains significantly triggered pine plants, specifically in primary metabolism that leads to growth promotion, which may help in survival rates after transplanting from the nursery to the field.In addition to nutrient mobilization, the auxin production trait should not be left out.Auxins are able to modify root structure, increasing root length and surface, thereby creating a higher surface to increase nutrient absorption capacity and fulfill the needs created on the shoot sinks [42,43].Although the root architecture was not studied, this fact has already been demonstrated by our group on other occasions, in which auxin-producing bacteria were able to modify the root surface significantly, which had a positive effect on the aerial growth of the plant [44][45][46].It is well known that many PGPR influence plant metabolism and can improve plant growth and fitness to better cope with different stress situations [9,47].
In addition to growth effects, the metagenomic analysis of the rhizosphere was performed in two periods: September 2019 (after six inoculations, autumn sampling time, before transplanting) and June 2020 (nine months after transplanting, spring sampling time) to explore whether any changes in biological diversity had taken place following inoculations and to establish a relationship with growth parameters.The Shannon index revealed no significant differences in the alfa biological diversity among the treatments (Figure 3a,b), indicating that inoculation does not significantly affect bacterial biological diversity.
Similarly, when analyzing the number of unique taxa, most of the taxa were present in the rhizospheres of all the treatments (Figures S1 and S2).The PCoA showed a greater grouping based on the treatments in the first sampling time (autumn, before transplanting) than in the second sampling time (spring, after transplanting).Taking into account the limitations of the analysis due to the low variance, as explained, and not wanting to reveal any seasonal effects, this seems to indicate that in autumn the inoculations generated differences in the microbial communities (beta diversity), but these differences were diluted with time and when transplanting the plants to the field.Despite this, inoculation with the different bacteria and the consortium had an effect on the rhizosphere microbial community since in all the cases (all the treatments and the control) there were exclusive taxa for each treatment, and the number of unique species was maintained at both sampling times, that is, it was maintained from pot to field (Figures S1 and S2).The rarefaction curves indicated that the maximum diversity had been reached; so, the presence of these unique taxa can only be due to the treatment.
The exponential increase in the use of biofertilizers in recent years is leading many scientists to question the effect of inoculation on ecosystem stability in general and on resident microbial communities prior to inoculation in particular [48].Omics techniques are making it possible to study this aspect in considerable depth, with both very marked impacts on the composition of resident communities and virtually no impact, as in our case [11,49].It is not clear from the studies carried out whether there is a relationship between the modifying capacity of the inoculations and the intended effect on the plants treated with the biofertilizers, but further research is needed to determine how this may affect the functioning of the biogeochemical cycles that could be altered by PGPR inoculations in the future [48].
In the current study, the predominant bacterial phylum at the two sample times was Proteobacteria, followed by Bacteroidetes, Acidobacteria, and Actinobacteria (Figure S3).This is not surprising, as these four groups are highly abundant in soil and play an important role in the recycling of nutrients through biogeochemical cycles.Proteobacteria have already been identified in other studies as the most prevalent bacterial group in the rhizosphere [50][51][52][53]; among other reasons, this is because of their decisive role in the nitrogen cycle [54].Bacteriodetes and Acidobacteria are two phyla of great importance in the recycling of organic matter and have an important role in the carbon cycle [55][56][57][58][59][60].Actinobacteria (Gram-positive soil-dwelling bacteria) are widely distributed and have a high degree of diversity [61].The production of secondary metabolites and the synthesis of cellulase and chitinase are the factors that contribute to their widespread distribution [62][63][64].
Interestingly, at lower taxonomic levels, such as class and order, we can clearly see the differential effect that each PGPR or the consortium has had on the four phyla described above, especially in the spring sampling, where the plants had to select bacteria existing in the natural soil.Z7.15 is the PGPR that produced the most relevant changes at the class and order level of the microbiome.For example, it increased the abundance of bacteria of the Alpha and Gamma Proteobacteria classes in its rhizosphere.Within the Alphaproteobacteria, it increased the abundance of bacteria of the order Rhizobiales, and within Gammaproteobacteria, it increased the bacteria of the order Pseudomonadales.The order Rhizobiales includes nitrogen-fixing bacteria, and the order Pseudomonadales contains a large number of bacteria described as PGPR [65].This treatment (Z7.15) also enhanced the abundance of bacteria of the order Actinomycetales, which can form symbiotic nitrogen-fixing associations with over 200 species of plants and can also serve as growthpromoting or biocontrol agents.Therefore, inoculation with Z715 bacteria may promote soil health by increasing nitrogen-fixing bacteria and PGPR bacteria.These and other examples can be observed in Figures S4-S6, demonstrating the decisive interaction that takes place between the inoculated PGPR and the plants, which are capable of determining the microbiome that inhabits their rhizosphere and which will play a decisive role in the survival of the plants and their acclimatization to different situations.
The linear discriminant analysis (LDA), which looks for significant differences in taxa composition between treatments, does indicate differences between treatments (Figure 4).In general, different taxonomic groups corresponding to Gram-negative bacteria were the most represented, except for the rhizosphere of the plants inoculated with Z4.3 before transplanting (autumn) where Gram-positive bacteria, especially Bacillus, formed the most represented group, which is not strange since Z4.3 is a Bacillus, although this relationship is no longer observed in spring.Within the Gram-negative bacteria, the Caulobacteraceae family (Z5.4-A) and the Sphingomonadales were abundant; the latter was especially abundant in several of the treatments (Z5.4-S,Z4.3-A, and CONS-S).The Sphingomonadales group is related to hydrocarbon-contaminated sites [66] so, they are bacteria with a wide use in the bioremediation of soils contaminated with this type of xenobiotics.Some of the significantly abundant taxa are related to serious opportunistic infections of humans, such as the genera Legionella (Z7.15 after transplanting) and Massilia (consortium before transplanting) and some species of the Sphingomonas genus (consortium in spring).The incidence of infections in humans brought on by opportunistic microorganisms has skyrocketed in recent years.The rhizosphere is one natural reservoir for opportunistic pathogens.This ecosystem is a "microbial hot-pot" because of its high nutritional content, which increases bacterial abundances, including those with strong antagonistic features [67].Numerous investigations have demonstrated that pathogenicity in humans and beneficial interactions with plants are caused by comparable, if not identical, activities.For instance, siderophore uptake systems or extracellular enzymes are common systems used by both human pathogens and beneficial bacteria.In the rhizosphere of plants inoculated with Z7.15 in spring, the genus Chthoniobacter was significantly represented.The genus Chthoniobacter is described in the literature as a genus capable of degrading many complex organic compounds, and it is therefore important in the proper functioning of the carbon cycle [68].
A key aspect of the study of microbial communities is the relationship between the structure and function of a microbial community within the ecosystem [69].Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) v1 is the first developed tool to predict potential functional genes from 16S rRNA metabarcoding and has been the most popular since its launch in 2013 [70].The method has been used to research different microbiomes [71][72][73][74], despite being verified using data from the Human Microbiome Project.
In this work, we used the pipeline EasyMAP, which uses PICRUSt to predict microbial functional content and LEfSe.The linear discriminant analysis (LDA) seeks significantly overrepresented functions among the treatments.The functions significantly associated with rhizosphere communities were influenced by the different treatments (Figure 5).Z7.15 in autumn was the treatment with the higher number of represented functions, including carbohydrate metabolism and secondary metabolite synthesis, among others.Functions related to neurodegenerative diseases, such as Parkinson's disease (Z5.4 after transplanting), Alzheimer's disease (Z7.15 after transplanting), or Huntington's disease (control after transplanting), also appeared in several treatments.There is increasing evidence of the effect of the microbiome and certain bacteria in it on the triggering or worsening of these types of diseases [75,76].It is not uncommon for such bacteria to occur in the rhizosphere, as it is well known that this ecosystem harbors extraordinary taxonomic and functional diversity [77,78].
In a recent review by Kong and Liu [11] on the effects of PGPR inoculation on the structure and function of rhizosphere microbial communities, the enormous importance that these changes can have on the effectiveness of PGPR is clearly described.For this reason, studies such as this one, in which a thorough study of the microbial communities after PGPR inoculation is carried out, is key to a good understanding of the plant growth promotion (PGP) mechanisms of PGPR.Understanding these modifications and trying to understand their impact is essential for the establishment of a strategic plant-microbe partnership to improve plant health and fitness in rapidly changing environments.

Conclusions
Based on the above, the selection capacity of plants in relation to the rhizobacteria that inhabit their roots has been verified once more, since the isolates were able to produce siderophores and solubilizing phosphates, activities that can help them to survive in the degraded soil in which pines live.
The use of PGPR proves once again to be a very powerful biotechnological tool in reforestation programs as the selected rhizobacteria clearly improved plant development, especially in compromised situations for plants, such as development in degraded soils.
Although neither alpha nor beta diversity was affected by PGPR inoculation, a deeper analysis of the microbial communities shows that there was an effect of inoculation as unique taxa were detected in the different treatments; taxa were significantly more represented in some treatments than in others; and functions were significantly more represented in the different treatments.In particular, inoculation with Z7.15 PGPR produced the most relevant changes at the class and order level of the microbiome, enhancing the abundance of bacteria belonging to the Rhizobiales and Pseudomonadales orders, which could promote soil health by increasing nitrogen-fixing bacteria and PGPR bacteria.

Figure 1 .
Figure 1.Percentage of strains with different PGPR capacities in each zone.

Figure 3 .
Figure 3. (a) Shannon diversity index (alfa diversity) of rhizosphere microbial community in the different treatments; (b) Shannon diversity index (alfa diversity) of rhizosphere microbial community in the both sampling times, before (nursery, autumn) and after transplanting (spring).

Figure 4 .
Figure 4. Graphics of linear discriminant analysis (LDA) effect size (LEfSe).Horizontal bars represent the effect size for each taxon.The length of the bar represents the log10 transformed LDA score, indicated by vertical dotted lines.The threshold on the logarithmic LDA score for discriminative features was set to 2.0.The taxon of bacteria with statistically significant change (p < 0.05) in the relative abundance is written alongside the horizontal lines.The name of the taxon level is abbreviated as p-phylum; c-class; o-order; f-family, and g-genus.Autumn sampling time: nursery, before transplanting.Spring sampling time: after transplanting.

Figure 4 .
Figure 4. Graphics of linear discriminant analysis (LDA) effect size (LEfSe).Horizontal bars represent the effect size for each taxon.The length of the bar represents the log10 transformed LDA score, indicated by vertical dotted lines.The threshold on the logarithmic LDA score for discriminative features was set to 2.0.The taxon of bacteria with statistically significant change (p < 0.05) in the relative abundance is written alongside the horizontal lines.The name of the taxon level is abbreviated as p-phylum; c-class; o-order; f-family, and g-genus.Autumn sampling time: nursery, before transplanting.Spring sampling time: after transplanting.

18 Figure 5 .
Figure 5. Predicted microbial functions.EasyMAP uses PICRUSt to predict microbial functional content and LEfSe to visualize the results.The functional categories are based on the metabolic pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG).The length of the bar represents the log10 transformed LDA score, indicated by vertical dotted lines.The threshold on the logarithmic LDA score for discriminative features was set to 2.0.The predicted microbial function with statistically significant change (p < 0.05) is written alongside the horizontal lines.Autumn sampling time: nursery, before transplanting.Spring sampling time: after transplanting.

Figure 5 .
Figure 5. Predicted microbial functions.EasyMAP uses PICRUSt to predict microbial functional content and LEfSe to visualize the results.The functional categories are based on the metabolic pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG).The length of the bar represents the log10 transformed LDA score, indicated by vertical dotted lines.The threshold on the logarithmic LDA score for discriminative features was set to 2.0.The predicted microbial function with statistically significant change (p < 0.05) is written alongside the horizontal lines.Autumn sampling time: nursery, before transplanting.Spring sampling time: after transplanting.
: Number of unique bacterial taxons (species) in the different treatments in Autumn sampling (plants in nursery, before transplanting); Figure S2: Number of unique bacterial taxon (species) in the different treatments in Spring sampling (after trasplanting); Figure S3: OTUS abundance at phylum level.Data obtained from analysis done by EasyMap; Figure S4: OTUS abundance at class level: (a) phylum Proteobacteria; (b) phylum Bacteoridetes; (c) phylum Acidobacteria and (d) phylum Actinobacteria.Data obtained from analysis done by EasyMap; Figure S5: OTUS abundance at order level: (a) class Alphaproteobacteria; (b) class Betaproteobacteria and (c) class Gammaproteobacteria.Data obtained from analysis done by EasyMap; Figure S6: OTUS abundance at order level: (a) phylum Bacteoridetes; (b) phylum Acidobacteria and (c) phylum Actinobacteria.Data obtained from analysis done by EasyMap.Author Contributions: Conceptualization, J.A.L., B.R.S., A.G.-V.and F.J.G.-M.; methodology, J.A.L. and A.G.-V.; formal analysis, J.A.L. and A.G.-V.; validation, J.A.L., B.R.S., A.G.-V.and F.J.G.-M.; writing-original draft preparation, J.A.L. and A.G.-V.; writing-review and editing, J.A.L., B.R.S., A.G.-V.and F.J.G.-M.All authors have read and agreed to the published version of the manuscript.Funding: The results presented in this manuscript were funded by the project IDI-20180364: Investigation of combined grafting and inoculation techniques aimed at the generation of Pinus pinea seedlings that are super-producers of pine nuts and at the recovery of degraded soils.Financing by FEDER funds through the Pluriregional Operational Program of Spain 2014-2020 and the CDTI (Center for

Table 1 .
Plant height (cm) and stem diameter (cm) in June 2019 (2 months after the first inoculation with a total of 2 inoculations) and June 2020 (1 year after in transplanted trees with a total of 6 inoculations).Letters a, b indicate significant differences (p < 0.05) in June 2019, and letters x, y, z indicate significant differences (p < 0.05) in June 2020.