Aridity modulates belowground bacterial community dynamics in olive tree.

Aridity negatively affects the diversity and abundance of edaphic microbial communities and their multiple ecosystem services, ultimately impacting vegetation productivity and biotic interactions. Investigation about how plant-associated microbial communities respond to increasing aridity is of particular importance, especially in light of the global climate change predictions. To assess the effect of aridity on plant associated bacterial communities, we investigated the diversity and co-occurrence of bacteria associated with the bulk soil and the root system of olive trees cultivated in orchards located in higher, middle and lower arid regions of Tunisia. The results indicated that the selective process mediated by the plant root system is amplified with the increment of aridity, defining distinct bacterial communities, dominated by aridity-winner and aridity-loser bacteria negatively and positively correlated with increasing annual rainfall, respectively. Aridity regulated also the co-occurrence interactions among bacteria by determining specific modules enriched with one of the two categories (aridity-winners or aridity-losers), which included bacteria with multiple PGP functions against aridity. Our findings provide new insights into the process of bacterial assembly and interactions with the host plant in response to aridity, contributing to understand how the increasing aridity predicted by climate changes may affect the resilience of the plant holobiont.

Introduction

Scarcity and unpredictability of precipitations in drylands impose strong constraints on plants' physiology, affecting their biodiversity and ecosystem multifunctionality (Maestre et al., 2012; Valencia et al., 2015; Delgado‐Baquerizo et al., 2016; Nunes et al., 2017). The persistence of vegetation in such ecosystems is highly dependent on the capacity of the plant to use water, to limit evapotranspiration and to overcome drought (Nunes et al., 2017; Dörken et al., 2020; Welles and Funk, 2020; Zogas et al., 2020), as well as to establish beneficial interactions with the surrounding edaphic microbial communities (Vandenkoornhuyse et al., 2015; Trivedi et al., 2020). The interactions of microorganisms and plants start in the rhizosphere, where a subset of the edaphic microorganisms is selected (Hartmann et al., 2009; Schlaeppi and Bulgarelli, 2015; Van Der Heijden and Schlaeppi, 2015; Sánchez‐Cañizares et al., 2017); these microbes offer to the plant multiple functions, such as nutrient recycling and solubilization, litter decomposition and plant growth promoting (PGP) services (Caruso et al., 2011; Mapelli et al., 2012; Hassani et al., 2018; Trivedi et al., 2020). Once associated to the plant, the microorganisms act as an extension of the host genotype able to rapidly respond to environmental changes (Zilber‐Rosenberg and Rosenberg, 2008; Barnard et al., 2013; Hardoim et al., 2015; Rosenberg et al., 2016; Leung et al., 2020). For example, it has been demonstrated that the adaptive responses of plant to drought is mainly governed by the belowground microbial communities and the rapid adaptation of microbial community to stressor (such as drought) increased the fitness of the plant to the same stress (Lau and Lennon, 2012). The ability of the microbiome to buffer plant response to stressful environmental conditions has been observed in several ecosystems, revealing how the interaction between these two components (i.e., plant and associated microorganisms) consistently leads to better performance of the holobiont (Rodriguez et al., 2008; Marasco et al., 2012, 2016; Cherif et al., 2015; Vandenkoornhuyse et al., 2015; Bang et al., 2018; Vigani et al., 2019; Alsharif et al., 2020; Trivedi et al., 2020). This is also confirmed by the fact that bacteria and fungi that carry biofertilization, biopromotion and bioprotection capacities are consistently recruited by the host, defining a ‘core stress microbiome’ (Timm et al., 2018; Marasco et al., 2018a; Mosqueira et al., 2019).

Aridity is commonly thought and has been demonstrated to negatively impact such positive associations among plant and edaphic microorganisms by reducing the abundance and diversity of key microbial taxa typically associated with soil fertility, with the consequence to threat the ecosystem resilience and functioning (Maestre et al., 2015; Delgado‐Baquerizo et al., 2016; Neilson et al., 2017; Karray et al., 2020). As recently shown by Berdugo and colleagues (2020), these changes occur sequentially, starting with vegetation decline and a soil disruption phase, followed by the final systemic breakdown of the ecosystem. Such land degradation process could potentially lead to extensive microbial phylotype changes and local extinctions, drastically affecting the microbial networking and the functioning of terrestrial ecosystems (Delgado‐Baquerizo et al., 2020a). Considering that by 2100 more than 20% of the land will be progressively exposed to an increase in aridity expected to negatively affect both ecosystem multifunctionality (Delgado‐Baquerizo et al., 2016; Berdugo et al., 2020; Delgado‐Baquerizo et al., 2020b) and food production (Busby et al., 2017; Toju et al., 2018; Alsharif et al., 2020; Singh et al., 2020), the ecological understanding of microbiome assembly in the plant root system under aridity stress is pivotal to elucidate microbial recruitment patterns that may enhance holobiont resilience to desertification (Cheng et al., 2019; Rudgers et al., 2020).

In this study, we aimed to evaluate the recruitment and selection of bacterial microbiome in the plant root system under conditions of increasing aridity. We selected as a model plant the olive tree (Olea europea L.), a crop able to grow in a wide aridity range, from the semiarid ecosystems of the Mediterranean basin, up to the pre‐desertic regions of the North Africa (Tunisia, Israel, Egypt; faostat.com) and in the Arabic Peninsula (Al‐Ruqaie et al., 2016). We selected five locations devoted to olive tree cultivation in Tunisia along a latitudinal gradient of aridity encompassing higher, middle and lower arid bioclimatic zones (Ben Ahmed et al., 2007; Gargouri et al., 2012; Verner et al., 2018). We hypothesized that (i) aridity acts as an environmental filter, affecting the overall soil physical and chemical properties, and consequently, the biodiversity and functional potential of the bacterial community along the gradient; (ii) due to the paucity of nutrient and microbial resources, the soil reconditioning effect mediated by the olive root system should be amplified by the increase of aridity; (iii) given the effect of aridity on microbial communities' diversity, the microbiome assembly in the olive tree root system should be particularly affected by taxa adapted to aridity under the ongoing global changes.

Results

Aridity and chemical gradient along Tunisian olive regions

We analysed five regions devoted to olive tree cultivation over an almost 400 km North–South aridity transect that encompassed higher (Zaghouan), middle (Chraitia and Gafsa) and lower (Neffatia and Matmata) arid regions, characterized by different amounts of rainfall per year (Table 1; Figs S1 and S2). Olive tree cultivation in these regions are managed by the application of traditional ‘desert farming’ techniques, including low tree density, virtuous use of water for irrigation and application of organic fertilizers (Ben Abdallah et al., 2020) that allow the cultivation of plants also in very harsh conditions (Fig. S2). According to the physico‐chemical analyses performed on bulk soil and root surrounding soil (Table S1), the studied olive orchards were characterized by significantly different edaphic conditions (PERMANOVA: F 3,8 = 144.3, p = 0.001, Fig. 1A; Table S2A), that were consistently modified by the activity of the olive root system (F 4,16 = 14.51, p = 0.001; Table S2B). As expected, the olive rhizosphere presented respect to the bulk soil: (i) significantly higher concentration of nutrients (N, P, K and Mg), organic matter (OM) and carbon, (ii) significantly higher microbial activity and (iii) significantly lower pH (Table S3). The rhizosphere soil recondition process varied along the aridity gradient (F 3,16 = 25.8, p = 0.001; Fig. 1A). Notably, physico‐chemical variations in both bulk and root surrounding soils followed a positive relationship with the annual rainfall (bulk soil: R 2 = 0.69, p < 0.0001 and root surrounding soils: R 2 = 0.72, p < 0.0001; Fig. 1B and C), showing an effect of aridity on the soil. The rate of rainfall was positively correlated with soil OM, total N, total organic C, exchange sodium percentage, cation exchange capacity and soil respiration (Fig. 1D), all components that play key roles in the formation and persistence of soil fertility and water retention (Schlesinger and Andrews, 2000; Darilek et al., 2009).

Table 1

List of stations identified along the North–South aridity transect; locations (or nearest by city), GPS position, bioclimatic zones and aridity index are indicated.

Station	Location	Lat (N)	Lon (E)	Bioclimatic zones	Aridity index
A	Zaghouan	36°24′36.0″	10°09′36.0″	Higher‐semiarid	0.17
B	Chraitia	35°13′48.0″	10°02′24.0″	Higher‐arid	0.43
C	Gafsa	34°24′12.2″	8°41′34.8″	Higher‐arid	0.58
D	Matmata	33°33′59.2″	9°42′59.9″	Middle‐arid	0.71
E	Neffatia	33°15′32.9″	10°50′17.3″	Lower‐arid	0.75

Fig 1

Edaphic physico‐chemical conditions along the aridity transect.

A. Principal coordinates analysis (PCoA) visualizes the similarity among orchards' bulk soil and olive root surrounding soil along the aridity transect based on Euclidean distance matrix; note, soils from Matmata orchard were not available for physico‐chemical analysis.

B and C. Precipitation decay patterns of bulk soil and root surrounding soils, respectively, showing the relationships between differences in annual rainfall and the diversity of physico‐chemical property as Euclidean distance; the relationship is tested by a linear regression (R 2) with a significance probability estimate (p‐value).

D. Pearson correlation among rainfall and physico‐chemical components of soils; significant correlation is indicated with stars: *, p < 0.05; **, p < 0.01; ***, p < 0.001; RSS: root surrounding soil, CEC: cation exchange capacity, ESP: exchange sodium percentage, TN: total N, TOC: total organic carbon, OM: organic matter, DOC: dissolved organic carbon, SoilResp1: soil respiration expressed as mg CO2 per g of dry soil, SoilResp2: soil respiration expressed as mg CO2 per g of TOC.

Distribution of bacterial members associated with olive trees' root system

Out of the 2849 bacterial operational taxonomic units (OTUs) detected (Table S4), a limited number was classified as abundant (>1% of relative abundance: 1.7%, 0.5%, 0.2% and 0.3% of OTUs in root tissues, rhizosphere, root surrounding soil and bulk soil, respectively), while a grand majority of OTUs were rare (<0.1% of relative abundance: 86%, 93%, 92% and 91% of OTUs, respectively; Fig. S3A–D). We plotted occupancy of OTUs versus their average relative abundance and a positive relationship was detected in all compartments (black circles in Fig. 2A; relative abundance of single OTUs in Fig. S3E,F): while OTUs with low abundance had low occupancy, those with high abundance showed high occupancy and were detected across the entire aridity gradient, suggesting a common mechanism to shape the bacterial community in olive root systems and orchards. Nevertheless, we found that the OTUs showed a different occupancy–frequency distribution in the four compartments: root tissues had a bimodal pattern distribution strongly skewed to the left of the plot (i.e., most species found in a single sample), while the rhizosphere and root surrounding soils had an increased number of taxa that occurred in all the sites (Fig. 2A). Notably, bulk soils were the only ones showing specific spikes of occupancy in concomitance with location changes (Fig. 2A), indicating a strong local‐specificity of edaphic bacterial community hosted by olive orchards.

Fig 2

Occupancy–frequency distribution of OTUs in olive orchards across the aridity transect.

A. Number of OTUs is reported in function of their occupancy expressed as number of samples where they are detected; number of OTUs are indicated by the left y‐axis. Black circles indicate relative abundance of OTUs (percentage of total sequences) for each occupancy; relative abundances are indicated by the right y‐axis. The number of samples is given by the x‐axis; maximum occupancy, n = 15.

B. Bar charts show the relative abundance of the main bacterial classes associated with root tissues, rhizosphere, root surrounding soil of olive trees and orchard bulk soils along the aridity transect; bacterial taxa belong to classes with relative abundance <1% are reported as ‘Others’.

Recruitment and selection of olive tree root system associated bacterial communities from the soils

Alpha diversity (Shannon index and richness) significantly varied along the compartments (ANOVA, Shannon index: F 3,56 = 109.5, p < 0.0001 and richness: F 3,56 = 80.88, p < 0.0001; Table 2); a sequential decrement was observed from bulk soil to the rhizosphere and the soil surrounding the roots, except for the nutrient‐rich soil of Zaghouan subjected to a lower level of aridity, in which all the edaphic bacterial communities had similar alpha diversity (Table 2). The lowest alpha diversity values were consistently measured in the root tissues (Table 2) where the edaphic recruitment mediated by plant culminates (Van Der Heijden and Schlaeppi, 2015). The number of local‐OTUs shared among all four fractions ranged from 6% to 13%, and OTUs present only in the plant‐affected soil (rhizosphere and root surrounding soil) ranged between 39% and 59% of OTUs. All together these OTUs accounted for 84%–94% of the bacterial relative abundance (Fig. S4). Subsets of compartment‐specific OTUs were also found with a consistent decremental trend from root surrounding soil (specific‐OTUs, 5%–14%) to rhizosphere (2%–6%) and root tissues (0.6%–3%), suggesting that the root systems and the olive orchard soils act as specialized niches for specific taxa. Since the bacterial abundance in the rhizosphere was relatively higher than in the bulk soil (in average, 1.1 × 109 and 2.7 × 107 CFU g−1 of soil, respectively; Table S5), we believe that the microbial community assembly in olive tree soil‐endosphere continuum follows a selection‐amplification model (Wang et al., 2020). The niche partitioning mediated by the root system was clearly reflected in the taxonomic distribution of the major taxa (Fig. 2B). While phylogenetically heterogeneous communities were detected in the different edaphic compartments (rhizosphere, root surrounding soil and bulk soil), Actinobacteria, Alphaproteobacteria and Gammaproteobacteria dominated the endophytic communities.

Table 2

Alpha diversity of bacterial communities associated with olive tree root systems (root tissues, rhizosphere, root surrounding soil) and orchards' bulk soil.

Alpha diversity	Location	Root tissue	Rhizosphere	RSS	Bulk soil
Richness	Zaghouan	342 ± 51a	1398 ± 50b	1434 ± 11b	1390 ± 15b
	Chaiatia	269 ± 90a	1697 ± 116b	1920 ± 48bc	2057 ± 210c
	Gafsa	205 ± 86a	1225 ± 154b	1675 ± 344bc	1683 ± 124c
	Matmata	333 ± 23a	1383 ± 153b	1581 ± 431bc	1830 ± 17c
	Neffatia	213 ± 20a	877 ± 147b	1283 ± 416bc	1337 ± 156c
Shannon index	Zaghouan	3.32 ± 0.45a	5.92 ± 0.06b	5.94 ± 0.05b	5.93 ± 0b
	Chaiatia	2.94 ± 0.23a	5.82 ± 0.47b	6.6 ± 0.13c	6.51 ± 0.01c
	Gafsa	3.28 ± 0.5a	5.26 ± 0.25b	6.04 ± 0.47bc	6.11 ± 0c
	Matmata	4.45 ± 0.13a	5.68 ± 0.27b	6.24 ± 0.26c	6.47 ± 0.01c
	Neffatia	3.8 ± 0.23a	5.1 ± 0.1b	5.71 ± 0.45bc	5.7 ± 0.01c

Values are expressed as average ± standard deviation (n = 3). Different lowercase letters denote significant mean difference in alpha diversity among compartments based on the pairwise Tukey's test at p‐value  < 0.05.

Variation of olive tree's root system microbiome along the aridity gradient

We examined the similarity among the structure of bacterial communities associated with olive trees' root system and bulk soil (beta‐diversity), adopting Bray–Curtis (BC) similarity distance measure. The results revealed that the four compartments harboured significantly distinct bacterial communities, consistently detected in all the olive orchards (PERMANOVA: F 15,40 = 10.41, p = 0.001; Fig. 3A; Table S6) and explaining up to 47% of the variation observed. Despite the association between olive tree and bacteria was established in a similar manner in the five locations, the distribution pattern of samples in the ordination space of PCoA showed a differentiation along the aridity gradient with significant differences among locations (F 4,40 = 12.76, p = 0.001; variation explained, 26%; Table S6). In line with the high environmental heterogeneity of the aridity gradient—mainly associated with variation in precipitation—a significant decrease of similarity among bacterial communities was found as a function of rainfall decay within each compartment but with a different extent (ANCOVA: F 3,412 = 3.1, p = 0.028): the greatest rate of decay was observed in bulk soils, followed by root surrounding soil, rhizosphere and root tissues (Fig. S5).

Fig 3

Beta‐diversity and rhizosphere effect.

A. Principal coordinates analysis (PCoA) shows the diversity of bacterial communities according to compartment and location.

B. Heat map reports Bray–Curtis similarity values among bacterial community hosted by bulk soil and olive root system (RT: root tissues, RH: rhizosphere, RSS: root surrounding soil).

C. Depleted (left side) and enriched (right side) OTUs in the rhizosphere compared with bulk soil are reported for Neffatia (orange) and Zaghouan (green). Only OTUs which have a log2‐normalized average read count and p‐value <0.01 are included and plotted as dots. In the x‐axis are reported the log2‐fold change values, while in y‐axis the number of OTUs showing such log2‐fold change values.

D and E. Percentage of enriched and depleted OTUs (log2‐fold change value >2) in the rhizosphere and root surrounding soil (RSS) compared with bulk soil, respectively.

Olive tree's rhizosphere effect along the aridity gradient

We observed a significantly different level of the so‐called ‘rhizospheric effect’ along the aridity gradient: the rhizosphere bacterial communities from olive grown in Zaghouan—with higher annual rainfall and more fertile soil—were more similar to those of bulk soil (BC‐similarity, 68.8%), while the bacterial communities of Neffatia olive trees' rhizospheres and related bulk soils were more distant among each other's (BC‐similarity, 43.7%; Fig. 3B); Chraitia, Gafsa and Matmata showed an intermediate similarity among rhizospheres and bulk soils, with BC‐similarity of 63.1%, 63.8% and 53.8% respectively. Although the rhizosphere effect showed different levels of intensity along the aridity gradients, the separation between bulk soil and rhizosphere in the PCoA ordination space was mainly along the axes PC2 for all the olive trees (Fig. S7), confirming that plant recruitment and selection processes in the rhizosphere are established in a similar manner but with different extent.

We investigated if the variation in the magnitude of ‘rhizosphere effect’ could be attributed to a differential OTU's enrichment/depletion between bulk soil and rhizosphere. In line with the multivariate analysis, the number of bacterial OTUs significantly enriched/depleted in the Zaghouan rhizosphere was lower than in orchards with reduced rainfall, such as the one of Neffatia and the other locations, along the whole gradient of fold‐changes tested (Fig. 3C and D). For instance, at a threshold of log2‐fold change (adjusted p‐value < 0.001), we found that 13.7% and 9.2% of bacterial OTUs were affected (enriched and depleted) in Neffatia and Zaghouan olive rhizospheres respectively (Fig. 3D). This trend was also observed in the root surrounding soil with a moderate extent (Fig. 3E).

Aridity‐winners and aridity‐losers: how aridity affects the plant–soil feedbacks

We investigated the distribution of OTUs along the aridity gradient by evaluating their relative abundance in function of annual rainfall. We found that the distribution of 121 and 381 bacterial taxa was negatively (Spearman‐ρ < −0.5) and positively (ρ > 0.5) correlated with precipitation, respectively (Fig. 4A and B; Table S7), as the sum of their relative abundance (R 2 = 0.83, p < 0.0001 and R 2 = 0.66, p < 0.0001; Fig. 4C and D). These relations revealed that taxa colonizing olive trees' root systems and orchard bulk soils have different environmental preferences (i.e., amount of rain) along the aridity gradient, acting as aridity‐winners and aridity‐losers, respectively (sensu from Delgado‐Baquerizo et al., 2020a). In particular, among the aridity‐winners, we found multiple OTUs belonging to Actinobacteria (36%), Thermoleophilia (23%), Chloroflexi (11%), Alphaproteobacteria (8%) and Bacilli (5%) classes (Fig. 4A; Table S7). On the contrary, the aridity‐losers (i.e., poorly adapted to low precipitations) were mainly affiliated to Acidobacteria (20%), Phycisphaerae (16%), Rubrobacteria (11%), Alphaproteobacteria (9%) Blastocatellia (8%), Gemmatimonadetes (8%) Gammaproteobacteria and Verrucomicrobiae (5% each one; Fig. 4B; Table S7). Functional profiles associated with the members of the aridity‐winner and aridity‐loser categories were also inferred by the bacterial 16S rRNA gene sequence data. The variations of predicted metabolic pathways was significantly different among the two categories (manyglm: F 1,88 = 139.8, p = 0.001), with the xenobiotics (p = 0.001), other secondary metabolites (p = 0.001), terpenoids and polyketides (p = 0.003), and glycan (p = 0.007) metabolisms explaining such difference. Considering only PGP traits, such as enzyme‐encoding genes putatively involved in activities of biofertilization [nitrogen metabolism, phosphate solubilization and siderophore synthesis] and biostimulation [auxin production, 1‐aminocyclopropane‐1‐carboxylic acid deaminase (ACCd) activity and volatile organic compounds (VOCs) production], aridity‐winner and aridity‐loser taxa carried distinct traits (F 1,88 = 49.91, p = 0.001); biofertilization activities were significantly enriched among aridity‐loser bacteria (Mann–Whitney test, p < 0.05), while biostimulation related traits were equally distributed among the two bacterial categories (p > 0.05).

Fig 4

Aridity‐winner and aridity‐loser bacteria associated with olive trees' root systems and orchard soils.

A and B. Taxonomic affiliation of bacterial OTUs negative (Spearman‐ρ < −0.5, p < 0.01) and positive (ρ > 0.5, p < 0.01) correlated with rain, therefore indicated as aridity‐winners and aridity‐losers, respectively.

C and D. Relation between precipitation (rain, mm) and total relative abundance of aridity‐winners and aridity ‐loser, respectively.

E and F. PGP potential and abiotic stress resistance of cultivable bacteria isolated from Neffatia and Zaghouan respectively. Taxonomic affiliation and aridity category are also reported.

We further investigated the PGP potential of cultivable bacteria isolated from the two extreme sites of the aridity gradient, Zaghouan and Neffatia. We isolated a total of 264 and 208 bacterial strains from the olive tree root systems (i.e., root tissues, rhizosphere and root surrounding soil) and orchards bulk soil in Zaghouan and Neffatia respectively (Table S8). By comparing the 16S rRNA gene sequences of the bacterial isolates with those of OTUs, we detected that 32% and 50% of bacterial isolates can be categorized as aridity‐winner in Zaghouan and Neffatia, respectively, while the remaining portion was generalist; any of the bacterial isolates were categorized as aridity‐loser category (Fig. 4E and F). This result was in accordance with the strategy of cultivation used in which oligotrophic media favoured the selection of stress‐tolerant bacteria, such as Actinobacteria (Arthrobacter, Microbacterium and Streptomyces) and Bacilli (Bacillus, Paenibacillus and Terribacillus), while rich media allowed the proliferation of generalist bacteria, including Proteobacteria and Bacteroidetes. The PGP potential of the isolates has been also evaluated by in vitro tests, attributing a final PGP score based on the number and level of PGP‐related activities. Bacterial isolates showed a PGP‐score ranging from 4 to 15, with consistently higher PGP‐score values within the aridity‐winner category compared with generalist bacteria (Zaghouan: t 1,117 = 5.91, p < 0.0001 and Neffatia: t 1,113 = 3.13, p = 0.002). Aridity‐winners showed a wide range of PGP‐related functions, including release of auxin, solubilization of phosphate and iron, fixation of nitrogen, production of EPS and release of ammonia (Fig. 4E and F).

Effect of aridity on bacterial interactions within communities associated with the olive tree belowground

We built a single co‐occurrence network between the members of the bacterial communities associated with olive trees' root systems and orchard bulk soils along the entire aridity gradient. We found that bacteria were grouped into six main independent modules (Fig. 5A; Table S9). The six modules were composed by 50% bacterial nodes computed in the co‐occurrence network: 15.3% (Module F), 10.8% (Module D), 7.4% (Module A), 6.2% (Module C), 5.4% (Module E) and 5% (Module B). Based on the OTUs aridity distribution (aridity‐winners and aridity‐losers; Fig. 4A and B), we analysed the composition of bacterial nodes in each module. Module A was dominated by aridity‐winners (73%), Module B showed 50% of nodes categorized as aridity‐winners and 50% as aridity‐losers, while the remaining Modules C–F were mainly composed by aridity‐losers (72%, 93%, 95% and 100%, respectively; Fig. 5B). Such distribution was significantly related with the region annual rainfall: Modules A and B were negatively related with rain (aridity‐tolerant modules), while Modules D, E and F were positively related with it (aridity‐sensitive modules); Module C did not show any significant relation with precipitation (p = 0.13; Fig. S7). Within each module, we investigated the taxonomic diversity of members that may respond in a similar manner to the aridity (Barberán et al., 2012). For instance, aridity‐sensitive modules (D, E and F) were mainly composed by bacterial taxa belong to Acidobacteria (Acidobacteriia, Blastocatellia, Holophagae and Subgroup 6), Planctomycetes (Phycisphaerae and Planctomycetacia), Gemmatimonadetes (Gemmatimonadetes and Longimicrobia), Verrucomicrobiae, Bacteroidia, Alphaproteobacteria and Gammaproteobacteria (Fig. 5C), carrying important PGP activities relevant in soil‐fertility and plant fitness, such as P solubilization, N fixation and VOCs release (Fig. 5D), while within aridity‐tolerant modules (A and B) we mainly detect members of Actinobacteria (Frankiales and Micrococcales) and Thermoleophila (Solirubrobacterales and Gaiellales) involved in biopromotion activities [auxin production and VOCs release; Fig. 5C and D]. Based on these results, we supposed that any modification of the climatic conditions (e.g., drought extent) may affect the abundance of the aridity‐loser bacteria with possible disruption of their networking, i.e., aridity‐sensitive Modules D, E and F. However, it is important to note that within these modules (except F) we detected also a limited number of aridity‐winner bacteria that compose the rare microbiome of the less‐arid olive orchard of Zaghouan.

Fig 5

Co‐occurrence correlation network of the microbiome associated with olive trees and orchard ecosystems. A. Diagram of co‐occurrence bacterial network with nodes coloured by their belonging to the different modules. B. Characterization of the taxa within each module is also reported, along with (C) the node categorization in aridity‐categories (aridity‐winner and aridity‐loser) and (D) their predicted PGP services. Taxa are coloured following the colour‐code reported in Fig. 2.

Discussion

Soils of arid and semi‐arid regions—approximately one‐third of the planet's surface (Laity, 2009)—are characterized by water deficiency that restricts plant and microbial activity (Maestre et al., 2015; Neilson et al., 2017; Ochoa‐Hueso et al., 2018; Shi et al., 2020). Due to current changes in climate (i.e., reduced rainfall, increasing drought and temperature), these areas are projected to increase in size by the end of this century (Huang et al., 2015). These changes in water availability can have profound and drastic impacts on vegetation, soil microbiomes and ecosystem functioning (Berdugo et al., 2020). In this study, we evaluated the belowground plant–soil feedbacks—intended as interactions among plants, edaphic bacterial community and abiotic soil conditions—along a North–South aridity gradient in Tunisia using olive tree (O. europea L. cv. Chemlali) as model plant, selected for its tolerance to drought events and the large diffusion of its cultivation in arid and semi‐arid ecosystems throughout the Mediterranean basin and Arabian Peninsula (Gargouri et al., 2012). Here, we showed in a field‐based observational study that the bacterial communities associated with olive tree root systems and orchard bulk soils are significantly influenced by the level of aridity, offering new information regarding the process by which plant‐associated bacterial communities respond to rainfall changes in a desert‐farming landscape.

Arid and semi‐arid environments—receiving no more than 100–500 mm of annual rain—are characterized by low levels of soil nutrients and OM, poor soil structure, high salinity and reduced water‐holding capacity that result in frequent drought events and loss of fertility (Maestre et al., 2015; Naylor and Coleman‐Derr, 2018). Along the aridity gradient studied, the annual average precipitation received by the olive orchards was significantly related to several physico‐chemical soil attributes, such as pH, nutrient level, OM, cation exchange capacity and soil respiration; all of which are considered important attributes to assess soil health and fertility, especially when considering soil capacity to support plant growth (Mader, 2002; Bünemann et al., 2018). Along with the effect of aridity, orchard soils are consistently influenced and modified by the rhizosphere effect mediated by the olive trees (van Dam and Bouwmeester, 2016). This process modulates the plant–soil feedbacks by modifying the soil physico‐chemical conditions and creating a favourable niche for edaphic microorganisms (Hartmann et al., 2009). As expected, we found that the rhizosphere effect was an important driver of modifications in soil physico‐chemical conditions and microbial communities' diversity in the olive root rhizosphere and surrounding soil along the entire aridity gradient. However, the entity of the rhizosphere‐recondition was differentially displayed in the different transect sites sampled and was in particular amplified by the aridity: we detected that the rhizosphere microbiome of olive trees growing in the most arid region of Neffatia significantly differed from soil communities at a greater extent than those of olive trees in the semi‐arid region of Zaghouan. It is known that the rhizosphere effect has an important role in the plant response to environmental change (Gargallo‐Garriga et al., 2018; Hu et al., 2018; Zhalnina et al., 2018). For instance, during water limitation (e.g., 150 mm of annual rain in Neffatia), plants alter the structure of root system architecture and rhizosphere metabolomes (Preece and Peñuelas, 2016; Gargallo‐Garriga et al., 2018; Preece et al., 2018; Xu et al., 2018; Williams and de Vries, 2020). All these changes driven by water availability have complex effects on the microbial communities associated with soil and plants by altering their composition and diversity (Naylor et al., 2017; Santos‐Medellín et al., 2017; Fitzpatrick et al., 2018; Naylor and Coleman‐Derr, 2018). In particular, it was observed that drought and aridity drive the enrichment of Gram‐positive bacteria (oligotrophic phyla, such as members of Firmicutes, Chloroflexi and Actinobacteria) and the depletion of Gram‐negative (copiotrophic phyla, including Bacteroidetes, Planctomycetes, Acidobacteria and Verrucomicrobia) with a conserved pattern among plants species and soils from all continents (Maestre et al., 2015; Preece and Peñuelas, 2016; Neilson et al., 2017; Naylor and Coleman‐Derr, 2018; Karray et al., 2020). These shifts in the soil and plant root microbiomes (Figs 2 and 3) are mainly caused by changes in the taxa relative abundance, rather than their abolition/appearance (Naylor and Coleman‐Derr, 2018). Despite we used an amplicon sequencing approach and we are unable to show the absolute abundance of species occurrences, we found that the relative abundance of certain bacterial taxa is influenced by local precipitation: while 5% of taxa were negatively associated with precipitation along the aridity gradient (aridity‐winner taxa), 14% was positively related to it (aridity‐loser taxa). Based on their distribution and co‐occurrence pattern (Figs 4 and 5), we hypothesized that these two bacterial categories have different lifestyles and functions and are ad hoc selected based on the environmental conditions. For instance, aridity‐winners are dominated by Actinobacteria (Actinobacteria, Thermoleophilia and unclassified Actinobacteria), Chloroflexi and Firmicutes (Bacilli); they are known for their persistence in dryland, hyper‐arid desert and desert (Angel et al., 2013; Maestre et al., 2015; Schulze‐Makuch et al., 2018; Marasco et al., 2018b; Willing et al., 2020), and for their physiological capacity to survive under drought/aridity conditions; e.g., given by thicker cell wall, formation of endospores and synthesis of osmolytes (Raddadi et al., 2012; Stevenson and Hallsworth, 2014; Makhalanyane et al., 2015). On the contrary, among aridity‐losers we found taxa from multiple phyla, such as Acidobacteria (Acidimicrobiia, unclassified Acidobacteria, Blastocatellia), Planctomycetes (Phycisphaerae), Verrucomicrobiae, Gammaproteobacteria, along with Rubrobacteria and Gemmatimonadetes (Fig. 4). These bacteria are of fundamental importance to the quality and fertility of soil, such soil OM decomposition, degradation of plant‐derived polymers, nutrient cycling, plant biopromotion and protection against abiotic and biotic stresses (Holmes et al., 2000; Marasco et al., 2012; Rolli et al., 2015, 2017; Banerjee et al., 2016; Soussi et al., 2016; Dedysh and Ivanova, 2019; Ivanova et al., 2020). Notably, it was recently proposed that the structural variations in the soil and root‐associated microbial communities induced by aridity/drought are mainly aimed at the selection of microbial assemblages adapted to abiotic stress with the recruitment of beneficial microorganisms able to tolerate perturbation (Xu et al., 2018) and enhance host tolerance to water stress (Goh et al., 2013; Fitzpatrick et al., 2018; Araya et al., 2020; Simmons et al., 2020; Williams and de Vries, 2020; Willing et al., 2020). For example, Bacillus spp., Paenibacillus spp. and Actinobacteria selected in arid ecosystem are recognized to increase drought resistance in several plants, both arboreal (date palm, grapevine) and herbaceous species (pepper, soybean, tomato, maize), through production of plant hormones and biochemical activity that help to mitigate water stress (Marasco et al., 2012; Yaish et al., 2015; Soussi et al., 2016; Vigani et al., 2019; Alsharif et al., 2020; Ayangbenro and Babalola, 2021).

Taxa distribution along the aridity gradient highlights the risk that the plant holobiont can be negatively affected by the predicted scenario of climate changes, i.e., decrease in precipitation, increase in temperature and frequency of drought events (Huang et al., 2015), especially in those environments that are already semi‐arid/arid, such as the sites in Tunisia studied here (Verner et al., 2018), because taxa may change their abundance and distribution based on their aridity preferences. Following the distribution pattern described here, the members of aridity‐loser category are more sensitive to disturbances (i.e., decrease in precipitation) and will tend to reduce their abundance under increasing aridity. As well, aridity can change their relative contribution to the network by affecting aridity‐sensitive modules (Barberán et al., 2012); these repercussions on microbial network stability may further impact plant–soil feedbacks and plant fitness by affecting the provision of essential ecosystem functions and services (Kaisermann et al., 2017; Neilson et al., 2017; de Vries et al., 2018).

Despite it is possible suppose that ecosystems dominated by aridity‐loser taxa (Zaghouan) may be more prone to stressful conditions due to climatic events than other sites, such as Neffatia, it is important to underline that bacterial communities associated with olive tree root system are also characterized by a great portion of taxa with low relative abundance (‘rare biosphere’; Fig. 2), including certain aridity‐winner taxa. As observed in several ecosystems, the presence of a highly diverse rare biosphere is typical of resilient bacterial communities in which the ecosystem functioning redundancy is conferred by the many species having intermediate/low relative abundances, which can potentially overcome possible environmental changes and stressors that can be suboptimal or detrimental for the dominant species (Lynch and Neufeld, 2015; Jousset et al., 2017). These dynamics, in a condition of further increasing aridity, could induce the aridity‐winner rare taxa to be vicariant and replace the loss of dominant members (aridity‐losers), ultimately maintaining the provision of the necessary ecosystem functions to the plant holobiont under intensifying environmental changes and stresses. This suggests a possible form of resilience in the plant‐associated microbiome that represents an important resource for olive tree cultivation under exacerbated climatic conditions.

While it is pivotal the protection of the environment from desertification, the knowledge of selection dynamics of drought‐resistant bacterial communities and their interaction with plant can identify suitable nature‐based solutions to favour plant resilience and soil multifunctionality, especially in arid and semi‐arid ecosystems. Our study shows that aridity determines changes in the plant‐microbial community composition and interactions by modulating the distribution of aridity‐tolerant (winner) and aridity‐susceptible (loser) bacterial taxa. Taking into account the role of bacteria in the ecosystem functioning of terrestrial ecosystems (among others, carbon and nutrient cycling, and plant productivity), any modification in the community assembly and interactions due to increasing aridity can have drastic repercussion on the entire ecosystem.

Experimental procedures

Site description and sampling

A North–South aridity transect was identified in Tunisia. It was composed by five locations that encompass different climatic conditions (Table 1; Fig. S1); aridity values were obtained from the package R Envirem (Title and Bemmels, 2017) while annual precipitation from the WordClim dataset (Fick and Hijmans, 2017) both at the resolution of 1 km2 and mapped using the package raster ad ggmap (Kahle and Wickham, 2013) in R. Monthly temperature (°C) and rainfall (mm) were obtained from climate knowledgeportal.worldbank.org (average data from 1991 to 2016) and visualized using heat‐map in GraphPad Prism 8. Within each location, traditional olive orchards were identified (Fig. S2) and three olive trees (O. europea L. cv. Chemlali) were randomly selected as representative of the cultivation; the three plants shared the same soil type within the orchard. As indicated by the owner of the plantations, all the sampled plants were in production with age ranging from 20 to 30 years and were self‐rooted. The permission for samples' collection was obtained by the private owner of the olive orchards. Occurrences of olive tree cultivation in Tunisia were collected from GBIF database (https://www.gbif.org) and visualized in R using the package rgbif (Chamberlain et al., 2016).

Samples of olive roots were collected at 20–40 cm depth where the root system was denser. The soil that surrounded the root system and bulk soil 4 m far from olive trees were also collected. All soils and roots samples were collected under sterile conditions using sterile tools. Recovered samples were stored at −20°C for chemical and nucleic acid analyses, and at 4°C for microbiological isolation.

Soil physico‐chemical analysis

The olive trees' root surrounding soil and orchards' bulk soils were collected and bring in the laboratory for physico‐chemical analyses. Soils were air‐dried, sieved at 20 mm and then grinded to 2 mm, and further characterized for total organic carbon (TOC) and total nitrogen (TN) contents, available phosphorous (P), cation exchange capacity (CEC), exchange sodium percentage (ESP) and pH by using soil science procedures (Faithfull, 2002). Exchangeable cations (K, Ca, Mg, Na) were extracted by BaCl2‐triethanolamine solution (pH 8.1) and determined by inductively coupled plasma mass spectrometry (Varian, Fort Collins, USA). OM was calculated from TOC. Soil respiration was measured by trapping with alkali the CO2 produced during soils incubation at 20°C in the laboratory for 21 days (ISO 16072:2002 Soil quality—Laboratory methods for determination of microbial soil respiration). Final data were reported both on the basis of dry soil (soil respiration 1) and TOC (soil respiration 2). For dissolved organic carbon (DOC), 100 g of soil was extracted with 200 ml of deionized water (soil/water ratio of 1:2, w/v) at 20°C for 30 min under agitation (130 times min−1). After the extraction, samples were centrifuged for 15 min at 6500 r min−1. Then supernatants were filtered with 0.45 μm Millipore membrane (Advantec MFS, Pleasanton, CA) and DOC determined by dichromate method. All the analyses were performed in triplicate. Note, Matmata soils were not available to perform physico‐chemical analyses. The physico‐chemical tables containing the data from bulk soils and root surrounding soils of the four sites were normalized and used to create a resemblance matrix using the Euclidean distance in PRIMER (Anderson et al., 2008). We used such matrix to perform the Principal Coordinate Analysis (PCoA) and to perform correlation with rainfall decay. Similar percentage (SIMPER) analysis was run to identifying the contribution of the soil physico‐chemical variables to the Euclidean distance metric between root surrounding soil and orchard bulk soil. Correlation among physico‐chemical variables and annual rainfall was also evaluated (Person, p‐value <0.01).

DNA extraction

In the laboratory, the portion of root system collected was processed in order to separate the rhizospheric soil from the root tissues as described in Mapelli and colleagues (2018). Briefly, the roots were placed in sterile tubes with 9 ml of physiological solution (9 g l−1 NaCl), vortexed for 5 min and centrifuged at 4000 rpm for 5 min; the supernatant was discarded and the rhizosphere collected. A total of 0.5 ± 0.1 g of the rhizospheric soil, root surrounding soils and bulk soils were used to extract the total DNA, using the FastDNA SPIN Kit for Soil (BIO 101 Systems Q‐BIO gene; CA, USA) following the manufacturer instruction. In case of the root tissues, they were first sterilized with ethanol and hypochlorite (Cherif et al., 2015) and further smashed with liquid nitrogen. One (±0.2) g of root tissues powder was used to extract the DNA, adopting the hot CTAB procedure using 1 g of grind root sample (Khan et al., 2007). The quality and the size of the extracted DNA were checked by electrophoresis on 0.8% agarose gels. DNA was quantified using NanoDrop 1000 spectrophotometer (Thermo Scientific, Waltham, MA).

DNA amplification, library preparation and sequencing

Bacterial 16S rRNA gene fragments (~450 pb) were PCR‐amplified using primers 375F and 795R primers targeting the hypervariable V3‐V4 regions at Macrogen (South Korea). Raw sequences were analysed with QIIME pipeline, including paired‐end merging, quality filtering, trimming and dereplication, of the sequences (Caporaso et al., 2010) as described by Booth and colleagues (2019). Chloroplast and mitochondria‐classified OTUs were discarded and non‐prevalent OTUs with relative abundance <0.001% were filtered out. We obtained a total of 1 165 425 high‐quality sequences from the root system compartments (root tissues, 128 051; rhizosphere, 345 824 and root surrounding soil, 365 263) of olive trees and bulk soils (326 287; Table S4). These reads were clustered into a total of 2849 bacterial OTUs (cut‐off of 97% as sequence identity), with bulk soil, rhizosphere, root surrounding soil and root tissues harbouring 2654, 2807, 2815 and 1010 OTUs respectively (Table S4). Raw reads have been deposited in the Short Reads Archive of NCBI under BioProject ID PRJNA748566.

Bacterial diversity and statistical analysis

The bacterial OTUs table was used to calculate alpha diversity indices (Shannon diversity and richness) in R using the Phyloseq package. Occupancy–abundance curves were generated by calculating the number of samples in which a certain OTU was detected and its total relative abundance, along with the number of OTUs present at each level of occupancy. The OTUs shared among different compartments in each olive orchard were defined by a Venn‐diagram analysis in R. Beta‐diversity of bacterial communities was analysed using the compositional BC similarity matrix of the relative log‐transformed OTUs‐table in PRIMER 6 (Anderson et al., 2008). The BC‐similarity matrix was used to perform the unconstrained analysis of PCoA and permutational multivariate analyses of variance (PERMANOVA) to statistically test the impact of each experimental factors: ‘Compartment’ (four levels: root tissues, rhizosphere, root surrounding root and bulk soil) and ‘Location’ (five levels: Zaghouan, Chraitia, Gafsa, Matmata and Neffatia) defined as fixed. PERMANOVA pair‐wise tests were also conducted to evaluate the effect of ‘Location’ and ‘Compartment’ in each site. Average similarity (%) was used to quantify the ‘rhizosphere effect’ driven by olive root system along the aridity gradient. OTUs were tested for differential abundance (enrichment and depletion; log2‐fold change) between rhizosphere and bulk soil, and root surrounding soil and bulk soil was calculated using DEseq2 package in R (Love et al., 2014). We finally evaluated the effect of precipitation on the relative abundance of all OTUs using Spearman correlation (p < 0.01) in GraphPad Prism 8; using the definition proposed by Delgado‐Baquerizo and colleagues (2020a), OTUs positively correlated with rain (Spearman‐ρ > 0.5; negatively correlated with aridity) were defined ‘aridity‐losers’, the ones negatively correlated with rain (Spearman‐ρ < −0.5; positively correlated with aridity) were defined ‘aridity‐winners’, and the ones that did not show significant correlation (p > 0.01 or Spearman‐ρ > −0.5 and <0.5) were considered generalist (i.e., their abundance is not influenced by aridity).

Bacterial PGP potential associated with olive trees and orchard soils

The OTUs table was used to predict the functional potential of bacterial communities using PICRUSt2 (Douglas et al., 2020). It generated a relative abundance of KEGG orthology (KO) groups associated with each sample depending on matches between the representative OTUs sequences and the KEGG organisms. Among KEEG pathways, we analysed the predicted genes related to metabolic pathways, along with the ones involved in the PGP mechanisms (Marasco et al., 2018a). To select the significant variable that accounts for the variability between winner and loser we did use a multivariate generalized linear analysis using the package R manyglm (Wang et al., 2012). In addition, bacterial cultivation was also conducted; 1 g of sterile root tissues, rhizosphere, surrounding root soil and bulk soil samples were shaken with 9 ml of sterile physiological solution (9 g l−1 NaCl). The suspensions were serially diluted and plated in triplicate on solidified King's medium B (KB), R2A medium (Oxoid) and R2A with 500 mM NaCl. The colony‐forming units (CFUs) for gram were determined. Cultivable bacteria from Zaghouan and Neffatia were further selected for identification and PGP screening; for each location, approximately 20 strains per medium per fraction were randomly selected, purified, dereplicated and identified, following the procedures described in Marasco et al. (2012). The sequences of the partial 16S rRNA genes for isolates were deposited in the GeneBank database under the accession numbers MZ600270–MZ600503. To evaluate the presence and role of bacterial isolates in the total community, 16S rRNA sequences of the isolates were blasted against aridity‐winner, aridity‐loser and generalist OTU sequences, using as threshold per cent of identity above 97% and coverage over 95% (Marasco et al., 2018a). The obtained strains were further test in vitro for their PGP activities, including phosphate solubilization, siderophore release, nitrogen fixation, ammonia production, auxin (indole‐3‐acetic acid) production, protease activity and exopolysaccharide release, along with tolerance to abiotic stresses (temperature: 4°C, 42°C and 50°C; salinity: 5%, 8% and 10% NaCl; osmotic stress: 10% and 20% of Polyethylene glycol) as described by Marasco and colleagues (2012).

Co‐occurrence network of bacterial community

We built a correlation network between the bacterial phylotypes that showed an average relative abundance higher than 0.01% by calculating all pairwise Spearman correlation coefficients among these bacterial taxa in Conet (Faust and Raes, 2016). We kept exclusively positive correlations with Spearman's correlation coefficient ρ > 0.5 and p < 0.01 in order to provide information on microbial taxa that may respond robustly and similarly to environmental conditions (i.e., aridity), as previously described (Barberán et al., 2012; Barberán and Casamayor, 2014; Delgado‐Baquerizo et al., 2018). This led to a network with 1000 co‐occurrences and 406 nodes (over 1447). The co‐occurrence network was visualized with Gephi (Bastian et al., 2009) and default parameters were used to identify modules of soil taxa strongly interacting with each other; only modules accounting for at least for 5% of OTUs were considered in the further analyses, while the remaining ones (n = 58) composed by few bacteria (average, 3.5 nodes each module) were not considered (full list in Table S9). Sum of relative abundances of OTUs belong to each module was calculated and used to evaluate module‐OTUs distribution in function of the annual rain (mm).

Fig. S1. Olive tree cultivation and climatic conditions along the aridity transect. Maps show (A) olive tree distribution in terms of presence of olive tree cultivation and (B) aridity index expressed as Thornthwaite aridity index, which relates annual moisture deficit to annual potential evapotranspiration; the five sampling locations (Zaghouan, Chraitia, Gafsa, Matmata and Neffatia) are indicated in the maps with red stars. Monthly (C) rainfall (mm), (D) maximum and (E) minimum temperatures (°C) are also reported by using heat maps.

Fig. S2. Olive orchards along the north–south aridity transect. Satellites images of the five locations selected for the sampling (A, Zaghouan; B, Chraitia; C, Gafsa; D, Matmata; E, Neffatia); in each location, the private olive orchard sampled is indicated with a red line.

Fig. S3. (A–D) Species rank–abundance curves for each compartment (i.e. root tissues, rhizosphere, root surrounding soil and bulk respectively). All compartments were characterized by a log‐normal distribution, with a few dominant OTUs and a tail of rare OTUs. (E,F) Occupancy–frequency distribution of OTUs in the four compartments. The number of sites sampled along the aridity transect is given in the x‐axis (maximum occupancy, n = 15); relative abundance (sum of reads in the entire dataset) is reported for each OTU in the y‐axis.

Fig. S4. Venn diagrams visualize the bacterial OTUs distribution across olive plant compartments (root tissue, rhizosphere and root surrounding soil) and orchards' bulk soil in the five locations; values in the Venn diagrams indicate the number of OTUs, while the values outside the number of OTUs not present in that specific location.

Fig. S5. Rain decay analysis showing the trend of bacterial community similarity in (A) root tissues, (B) rhizosphere, (C) root surrounding soil and (D) bulk soil. In the x‐axis is indicated the difference among the precipitation (rain, mm) measured at each location. Regression lines: (A) root tissues, slope = −0.036; R 2 = 0.1, p = 0.001; n = 105, DBC‐similarity = 20.3%; (B) rhizospheres, slope = −0.047; R 2 = 0.19, p < 0.0001, n = 105, DBC‐similarity = 35.4%; (C) root surrounding soils, slope = −0.048, R 2 = 0.21, p < 0.0001, n = 105, DBC‐similarity = 35.6%; (D) bulk soil, slope = −0.083, R 2 = 0.22, p < 0.0001, n = 105, DBC‐similarity = 55.4%.

Fig. S6. Principal coordinates analyses (PCoA) visualize the similarity among olive plant compartments (root tissue, rhizosphere and root surrounding soil) and orchards' bulk soil in the five locations based on Bray–Curtis distance matrices.

Fig. S7. Relationships between rain and the relative abundance of nodes within the six modules. R 2 and p‐value of regressions are reported in each graph.

Click here for additional data file.

Table S1. Results of soil physico‐chemical analyses are reported for root surrounding soil and bulk soil along the aridity gradient; note that soils from Matmata are not available for physico‐chemical analyses. Values are expressed as mean ± standard deviation of three replicates. Different lowercase letters denote significant mean difference among locations based on the pairwise Tukey's test at p < 0.05.

Table S2. Permutational multivariate analysis of variance (PERMANOVA) pairwise comparison of soil physico‐chemical datasets for bulk soil and root surrounding soil (RSS) along the aridity gradient; note that Matmata is not included in the comparison.

Table S3. Similar percentage (SIMPER) analysis identifying the contribution percentage of each physico‐chemical parameter to the Euclidean distance between root surrounding soil (RSS) and bulk soil.

Table S4. Reads and number of OTUs for the four compartments (root tissues, rhizosphere, root surrounding soil and bulks soil) across the five locations along the aridity transect. Values are reported for the three replicates analysed.

Table S5. Colony‐forming unit (CFU) per g of soil obtained from King's B, R2A and R2A 5% NaCl isolation media. Values are expressed as mean and standard deviation of six individual plates (two serial dilutions for each medium).

Table S6. Results of PERMANOVA main and pairwise tests for location and compartment as fixed factors. RSS, root surrounding soil.

Table S7. List of OTUs correlated with rain. Significant (p < 0.01) positive and negative correlations are considered based on Spearman correlation value >0.5 and <−0.5 respectively.

Click here for additional data file.

Table S8. Taxonomic identification and in vitro PGP activity of bacterial isolates from olive tree root system and orchards bulk soil in Zaghouan and Neffatia.

Click here for additional data file.

Table S9. Co‐occurrence network modularity. List of modules detected by Gephi within the co‐occurrence network.

Click here for additional data file.